############################################################
#
#   Sovereign Debt Ratchet and Welfare Destruction
#   January 2023
#
#############################################################

#======================================================================
# Setup
#======================================================================

rm(list=ls())           # Clear the workspace
set.seed(908)

# load libraries that are required
packages_req <- c('AER','akima','aod','car','coda','data.table','deSolve', 'DTMCPack', 
                  'dplyr','emdbook','gdata','gplots','ggplot2','jrvFinance',
                  'MASS','Matrix','mnormt','mvtnorm','nloptr','orthopolynom','pastecs',
                  'polynom','plot3D','pryr','RColorBrewer','rgl','rootSolve','Rsolnp',
                  'sandwich','scales','stargazer','systemfit','texreg','xtable')

# install packages not yet installed
installed_packages <- packages_req %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages_req[!installed_packages])
}

# load all required packages
lapply(packages_req, require, character.only = TRUE)

#==================
# Directories setup
#==================

# set working directory to the directory where code is saved
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# specify main directory
main_dir <- getwd()

# create a "dump" directory in case this does not already exists
if (!file.exists("Dump")) dir.create(file.path(main_dir, "Dump"))
dump_dir <- file.path(main_dir, "Dump")

# specify the R file from which main functions are taken
source('Sovereign_Debt_Functions.R')

# construct output directories -- where all output results will be stored
todays_date   <- gsub("-", "", Sys.Date())

# create output directory in case such directory does not exist yet
if (!file.exists(todays_date)) dir.create(file.path(main_dir, todays_date))

# store output folder into variable "output_folder"
output_dir    <- file.path(main_dir, todays_date)
setwd(output_dir)

# plotting parameters 
cex_plot   <- 1.5
thickness  <- 2

# numerical parameters (for all calculations except for Aguiar-Amador) 
num_param_0             <- c()
num_param_0['n_x']      <- 2000
num_param_0['dt']       <- .25
num_param_0['max_iter'] <- 1000
num_param_0['tol']      <- .00001

#==============================
# Parameters -- Base Case Model
#==============================

param_0              <- c()
param_0['mu']        <- .02     # income growth rate
param_0['sigma']     <- .1      # income volatility
param_0['m']         <- 1/20    # bond maturity intensity
param_0['theta']     <- .5      # debt haircut upon default
param_0['alpha']     <- .96     # GDP drop upon default
param_0['r']         <- .05     # international risk free rates
param_0['delta']     <- param_0['r'] + .02  # impatience of government
param_0['hat_delta'] <- .5*param_0['r']+.5*param_0['delta']  # impatience of citizens
param_0['kappa']     <- param_0['r']                         # coupon rate
param_0['nu']        <- 0.4     # international risk prices
param_0['eta']       <- 0       # tax on issuance proceeds

# plotting parameters
plot_0              <- c()
plot_0['n']         <- 1000
plot_0['thick']     <- thickness
plot_0['font_size'] <- cex_plot
plot_0['axis_size'] <- cex_plot

#======================================
# Plot smooth no-commitment equilibrium
#======================================

# create output directory for pictures and calculations of the benchmark risk-neutral case
setwd(output_dir)
if (!file.exists('base_case')) dir.create(file.path(output_dir, 'base_case'))

# store output folder into variable "output_folder"
analysis_dir    <- file.path(output_dir, 'base_case')
setwd(analysis_dir)

# plot all functions of interest
plot_fn(param = param_0, plotting_param = plot_0)

#====================================================
# Plot smooth no-commitment equilibrium when g(0) = 0
#====================================================

# create output directory for pictures and calculations of the case g(0) = 0
setwd(output_dir)
if (!file.exists('base_case_2')) dir.create(file.path(output_dir, 'base_case_2'))

# store output folder into variable "output_folder"
analysis_dir    <- file.path(output_dir, 'base_case_2')
setwd(analysis_dir)

# parametrize such that xi<2, i.e. such that delta+m<sigma^2-(delta+mu)
param_0b          <- param_0
param_0b['nu']    <- 0
param_0b['sigma'] <- ceiling(1.2*sqrt((param_0b['delta']+2*param_0b['m']+param_0b['mu']))/.05)*.05

# plot all functions of interest
plot_fn(param = param_0b, plotting_param = plot_0)

#===========
# Simulation
#===========

# create output directory for pictures for simulation and illustrative plots
setwd(output_dir)
if (!file.exists('simulation')) dir.create(file.path(output_dir, 'simulation'))

# store output in appropriate folder
analysis_dir    <- file.path(output_dir, 'simulation')
setwd(analysis_dir)

# modify parameters so that impatience is not too big
param_sim      <- param_0

# parameters of simulation
sim_0          <- c()
sim_0['T']     <- 20     # length of time of simulation
sim_0['dt']    <- 1e-4   # time step
sim_0['Y0']    <- 1      # initial GDP level
sim_0['F0']    <- .25    # initial debt-to-income level

# perform simulation and generate plots
sim_fn(param_vec = param_sim, sim_param = sim_0, plotting_param = plot_0)

#==========================================
# Compute citizens' value -- for range of m
#==========================================

# create output directory for pictures for debt maturity sensitivity
setwd(output_dir)
if (!file.exists('maturity_sensitivity')) dir.create(file.path(output_dir, 'maturity_sensitivity'))

# store output in appropriate folder
analysis_dir    <- file.path(output_dir, 'maturity_sensitivity')
setwd(analysis_dir)

# assume citizen have same discount rate as international creditors
hat_delta       <- param_0['r']
n_mat           <- 50
res_df          <- array(NA, dim = c(n_mat,4))
colnames(res_df)<- c('maturity','welfare_loss','erg_v','erg_hat_v')
res_df[,'maturity'] <- seq(from = 2, to = 40, length.out = n_mat)

# compute welfare loss for range of maturities
param_2               <- param_0
param_2["hat_delta"]  <- hat_delta
for (i in (1:n_mat)){
  print(i)
  param_2["m"] <- 1/res_df[i,'maturity']
  
  # compute welfare loss using hat delta 1
  hat_v_out       <- no_commit_hat_value_fn(param_2, num_param_0)
  erg_dens        <- no_commit_ergodic_dist_fn(param_2)
  x_temp_bar      <- no_commit_cutoff_fn(param_2)
  x_temp          <- seq(from = 0, to = x_temp_bar, length.out = 1000)
  res_df[i,'welfare_loss'] <- hat_v_out[["hat_v0"]]*(param_2["hat_delta"]-param_2["mu"]) - 1
  res_df[i,'erg_v']        <- erg_dens$moments['v','mean']
  res_df[i,'erg_hat_v']    <- sum(erg_dens$density(x_temp)*hat_v_out$hat_v(x_temp))/sum(erg_dens$density(x_temp))
}

# plot hat v0 - autarky value as function of m
pdf("citizen_welfare_loss_vs_m.pdf")
plot(x = res_df[,'maturity'], y = 100*res_df[,'welfare_loss'], type='l', xlab = '1/m (years)', 
     ylab = 'welfare loss vs. autarky (%)',
     col="blue", lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n')
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot hat v0 - autarky value as function of m -- create .eps file
setEPS()
postscript("citizen_welfare_loss_vs_m.eps")
plot(x = res_df[,'maturity'], y = 100*res_df[,'welfare_loss'], type='l', xlab = '1/m (years)', 
     ylab = 'welfare loss vs. autarky (%)',
     col="blue", lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n')
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot ergodic welfare as a function of m
pdf("ergodic_welfare_vs_m.pdf")
par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
plot(x = res_df[,'maturity'], y = res_df[,'erg_v'], bty = 'n', type = 'l',axes = 'F', main = '',  
     xlab = '', ylab = '', col = 'blue', lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
axis(2, col="black", lwd = 1, cex.axis = cex_plot)
mtext(2,text = "", line = 2, cex = cex_plot)
par(new=T)
plot(x = res_df[,'maturity'], y = res_df[,'erg_hat_v'], axes=F, xlab="", ylab="", 
     type="l", main="", col='red', lty=2, lwd = thickness)
axis(1, col="black", lwd = 1, cex.axis = cex_plot)
mtext(1,text='1/m (years)',line=2, cex = cex_plot)
axis(4, col="black", lwd = 1, cex.axis = cex_plot)
mtext(4,text = "", line = 2, cex = cex_plot)
legend("bottomright", legend = c(expression(E~"["~v(x)~"] (l.h.s.)"),
                                 expression(E~"["~hat(v)(x)~"] (r.h.s.)")), 
       col = c('blue', 'red'), bty="n", lty = c(1,2), cex = cex_plot, lwd = thickness)
dev.off()

# generate csv file
write.table(res_df, file = paste("results.csv",sep=""), eol = "\n", na = "NA", fileEncoding = "", row.names = F)

#===========================================================
# Compute citizens' value function -- for range of hat delta
#===========================================================

# create output directory for pictures for hat delta sensitivity
setwd(output_dir)
if (!file.exists('hat_delta_sensitivity')) dir.create(file.path(output_dir, 'hat_delta_sensitivity'))

# store output in appropriate folder
analysis_dir    <- file.path(output_dir, 'hat_delta_sensitivity')
setwd(analysis_dir)

# set up parameters
n_delta              <- 50
output               <- array(NA, dim = c(n_delta,4))
colnames(output)     <- c('hat_delta','hat_v0', 'n_iter', 'error')
output[,'hat_delta'] <- seq(from = param_0['r'], to = param_0['delta'] - 1E-4, length.out = n_delta)

# compute welfare loss for range of citizen's discount rates
param_3               <- param_0
for (i in (1:n_delta)){
  print(i)
  param_3['hat_delta'] <- output[i,'hat_delta']
  hat_v_out            <- no_commit_hat_value_fn(param_3, num_param_0)
  output[i,'hat_v0']   <- hat_v_out[["hat_v0"]]*(param_3['hat_delta']-param_3['mu']) - 1
  output[i,'n_iter']   <- hat_v_out[["n_iter"]] - 1
  output[i,'error']    <- hat_v_out[["error"]]
}

# fill the case delta = hat_delta
output                        <- rbind(output,rep(NA,4))
output[n_delta+1,'hat_delta'] <- param_3['delta']
output[n_delta+1,'hat_v0']    <- 0

# plot hat v0 - autarky value as function of hat_delta
pdf("citizen_welfare_loss_vs_hat_delta.pdf")
plot(x = output[,'hat_delta'], y = 100*output[,'hat_v0'], type='l', xlab = expression(hat(delta)), ylab = 'welfare loss vs. autarky (%)',
     col = "blue", lwd = thickness,  cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n',
     ylim = c(min(100*output[,'hat_v0']),0))
abline(h = 0, lty = 2, col = 'black', lwd = thickness)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot hat v0 - autarky value as function of hat_delta -- create .eps file
setEPS()
postscript("citizen_welfare_loss_vs_hat_delta.eps")
plot(x = output[,'hat_delta'], y = 100*output[,'hat_v0'], type='l', xlab = expression(hat(delta)), ylab = 'welfare loss vs. autarky (%)',
     col = "blue", lwd = thickness,  cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n',
     ylim = c(min(100*output[,'hat_v0']),0))
abline(h = 0, lty = 2, col = 'black', lwd = thickness)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# generate csv file
write.table(output, file = paste("results.csv",sep=""), eol = "\n", na = "NA", fileEncoding = "", row.names = F)

#================================================
# Sensitivity analysis
# Ergodic default rate and ergodic debt-to-income
#================================================

# create output directory for pictures for comparative statics
setwd(output_dir)
if (!file.exists('comparative_statics')) dir.create(file.path(output_dir, 'comparative_statics'))

# store output in appropriate folder
analysis_dir    <- file.path(output_dir, 'comparative_statics')
setwd(analysis_dir)

# min parameters for sensitivity analysis
n_param    <- 10
comp_param <- list()
comp_param[['delta']] <- seq(from = param_0['r']+.01, to = param_0['r']+.05, length.out = n_param)
comp_param[['mu']]    <- seq(from = param_0['mu']-.01, to = param_0['mu']+.01, length.out = n_param)
comp_param[['sigma']] <- seq(from = param_0['sigma'], to = param_0['sigma']+.1, length.out = n_param)
comp_param[['kappa']] <- seq(from = param_0['kappa']-.02, to = param_0['kappa']+.02, length.out = n_param)
comp_param[['theta']] <- seq(from = param_0['theta']-.1, to = param_0['theta']+.1, length.out = n_param)
comp_param[['alpha']] <- seq(from = param_0['alpha']-.02, to = param_0['alpha']+.02, length.out = n_param)
comp_param[['m']]     <- seq(from = 1/40, to = 1/5, length.out = n_param)
comp_param[['r']]     <- seq(from = param_0['delta']-.05, to = param_0['delta']-.01, length.out = n_param)
comp_param[['nu']]    <- seq(from = 0, to = 1, length.out = n_param)
num_sensi             <- length(comp_param)
sensi_list            <- as.list(rep(NA, num_sensi))
names(sensi_list)     <- names(comp_param)

# outcome to store -- as a function of state variable x
deft_rate         <- as.list(rep(NA, num_sensi))   # list of ergodic default rates
att_pt            <- as.list(rep(NA, num_sensi))   # list of attraction points
zero_drift_pt     <- as.list(rep(NA, num_sensi))   # list of zero drift points
zero_ca_pt        <- as.list(rep(NA, num_sensi))   # list of zero current account points
def_bdry          <- as.list(rep(NA, num_sensi))   # list of default boundaries

# outcome to store -- as a function of state variable z = (kappa+m)x
z_att_pt            <- as.list(rep(NA, num_sensi))   # list of attraction points
z_zero_drift_pt     <- as.list(rep(NA, num_sensi))   # list of zero drift points
z_zero_ca_pt        <- as.list(rep(NA, num_sensi))   # list of zero current account points
z_def_bdry          <- as.list(rep(NA, num_sensi))   # list of default boundaries

# declare moments
moments           <- as.list(rep(NA, num_sensi))

# names
names(deft_rate)  <- names(comp_param)
names(att_pt)     <- names(comp_param)
names(def_bdry)   <- names(comp_param)
names(z_att_pt)   <- names(comp_param)
names(z_def_bdry) <- names(comp_param)
names(moments)    <- names(comp_param)

# compute all equilibria for different parameter tests
for (param_type in names(comp_param)){
  param_sensi                 <- param_0
  deft_rate[[param_type]]     <- array(NA,n_param)
  
  # outcome in "x" state variable
  att_pt[[param_type]]        <- array(NA,n_param)
  zero_drift_pt[[param_type]] <- array(NA,n_param)
  zero_ca_pt[[param_type]]    <- array(NA,n_param)
  def_bdry[[param_type]]      <- array(NA,n_param)
  # outcome in "z = (kappa+m)x" state variable
  z_att_pt[[param_type]]        <- array(NA,n_param)
  z_zero_drift_pt[[param_type]] <- array(NA,n_param)
  z_zero_ca_pt[[param_type]]    <- array(NA,n_param)
  z_def_bdry[[param_type]]      <- array(NA,n_param)
  # moments
  moments[[param_type]]       <- array(NA, dim = c(n_param,9,4), 
                                       dimnames = list(c(1:n_param),
                                                       c('x','z','spread','cons','ca','cons_vol','g','v','d'),
                                                       c('mean','stdev','skewness','kurtosis')))
  
  for (i in 1:n_param){
    print(paste("parameter: ",param_type," - completion: ", round(100*i/n_param,0),"%",sep=""))
    param_sensi[param_type]        <- comp_param[[param_type]][i]
    dummy                          <- no_commit_ergodic_dist_fn(param_sensi)
    deft_rate[[param_type]][i]     <- dummy$default_rate
    # outcomes with "x" as state variable
    att_pt[[param_type]][i]        <- no_commit_attraction_pt_fn(param_sensi)
    zero_drift_pt[[param_type]][i] <- no_commit_zero_drift_pt_fn(param_sensi)
    zero_ca_pt[[param_type]][i]    <- no_commit_zero_ca_pt_fn(param_sensi)
    def_bdry[[param_type]][i]      <- no_commit_cutoff_fn(param_sensi)
    # outcomes with "z" as state variable
    z_att_pt[[param_type]][i]        <- att_pt[[param_type]][i]*(param_sensi['kappa']+param_sensi['m'])
    z_zero_drift_pt[[param_type]][i] <- zero_drift_pt[[param_type]][i]*(param_sensi['kappa']+param_sensi['m'])
    z_zero_ca_pt[[param_type]][i]    <- zero_ca_pt[[param_type]][i]*(param_sensi['kappa']+param_sensi['m'])
    z_def_bdry[[param_type]][i]      <- def_bdry[[param_type]][i]*(param_sensi['kappa']+param_sensi['m'])
    # moments
    moments[[param_type]][i,,]     <- array(unlist(dummy$moments), dim = c(9,4))
  }
  
  # generate csv file
  setwd(analysis_dir)
  output <- data.frame('parameter' = comp_param[[param_type]],
                       'default_rate' = deft_rate[[param_type]], 
                       'attraction_point' = att_pt[[param_type]],
                       'zero_drift_point' = zero_drift_pt[[param_type]],
                       'zero_ca_point' = zero_ca_pt[[param_type]],
                       'default_boundary' = def_bdry[[param_type]],
                       'z_attraction_point' = z_att_pt[[param_type]],
                       'z_zero_drift_point' = z_zero_drift_pt[[param_type]],
                       'z_zero_ca_point' = z_zero_ca_pt[[param_type]],
                       'z_default_boundary' = z_def_bdry[[param_type]])
  write.table(output, file = paste0("output_sensi_",param_type,".csv"), eol = "\n", na = "NA", fileEncoding = "", row.names = F)
  write.table(moments[[param_type]][,,'mean'], file = paste0("mean_sensi_",param_type,".csv"), eol = "\n", na = "NA", fileEncoding = "", row.names = F)
  write.table(moments[[param_type]][,,'stdev'], file = paste0("stdev_sensi_",param_type,".csv"), eol = "\n", na = "NA", fileEncoding = "", row.names = F)
}

# plot all relevant sensitivities for all parameter tests
setwd(analysis_dir)
for (param_type in names(comp_param)){
  param_sensi <- param_0
  
  # plot default rate and credit spread sensitivity
  file_name <- paste("deft_rate_and_spreads_sensi_",param_type,".pdf",sep="")
  y_max     <- 125*max(deft_rate[[param_type]], moments[[param_type]][,'spread','mean'])
  y_min     <- 100*min(deft_rate[[param_type]], moments[[param_type]][,'spread','mean'])
  pdf(file_name)
  par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
  plot(x = comp_param[[param_type]], y = 100*deft_rate[[param_type]], lty = 2, 
       bty = 'n', type = 'l',axes = 'F', main = '', xlab = '', ylab = '',col = 'blue', lwd = thickness, cex.lab = cex_plot, 
       cex.axis = cex_plot, ylim = c(y_min,y_max))
  lines(x = comp_param[[param_type]], y = 100*moments[[param_type]][,'spread','mean'], col = 'red', lwd = thickness, lty = 3)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  axis(2, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(2,text = "ergodic average (% p.a.)", line = 2, cex = cex_plot)
  axis(1, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(1,text=param_type, line=2, cex = cex_plot)
  legend("topright", legend = c('default rate', 'credit spread'), 
         col = c('blue', 'red'), bty="n", lty = c(2,3), cex = cex_plot, lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # plot debt-to-income sensitivity with "x" as state variable
  file_name <- paste("debt_to_income_sensi_",param_type,".pdf",sep="")
  y_max     <- 100*max(att_pt[[param_type]], zero_drift_pt[[param_type]], zero_ca_pt[[param_type]], def_bdry[[param_type]], moments[[param_type]][,'x','mean'])
  y_min     <- 100*min(att_pt[[param_type]], zero_drift_pt[[param_type]], zero_ca_pt[[param_type]], def_bdry[[param_type]], moments[[param_type]][,'x','mean'])
  pdf(file_name)
  par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
  plot(x = comp_param[[param_type]], y = 100*att_pt[[param_type]], lty = 1, 
       bty = 'n', type = 'l',axes = 'F', main = '', xlab = '', ylab = '',col = 'blue', lwd = thickness, cex.lab = cex_plot, 
       cex.axis = cex_plot, ylim = c(y_min,y_max))
  lines(x = comp_param[[param_type]], y = 100*zero_drift_pt[[param_type]], col = 'pink', lwd = thickness, lty = 2)
  lines(x = comp_param[[param_type]], y = 100*zero_ca_pt[[param_type]], col = 'orange', lwd = thickness, lty = 3)
  lines(x = comp_param[[param_type]], y = 100*moments[[param_type]][,'x','mean'], col = 'green', lwd = thickness, lty = 4)
  lines(x = comp_param[[param_type]], y = 100*def_bdry[[param_type]], col = 'red', lwd = thickness, lty = 5)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  axis(2, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(2,text = "debt-to-GDP (%)", line = 2, cex = cex_plot)
  axis(1, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(1,text=param_type, line=2, cex = cex_plot)
  legend("topright", legend = c('attract. point','zero drift','zero current act.','ergodic mean','default point'), 
         col = c('blue', 'pink', 'orange','green','red'), bty="n", lty = 1:5, cex = cex_plot, lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # plot debt-to-income sensitivity with "z" as state variable
  file_name <- paste("debt_coverage_ratio_sensi_",param_type,".pdf",sep="")
  y_max     <- 100*max(z_att_pt[[param_type]], z_zero_drift_pt[[param_type]], z_zero_ca_pt[[param_type]], z_def_bdry[[param_type]], moments[[param_type]][,'z','mean'])
  y_min     <- 100*min(z_att_pt[[param_type]], z_zero_drift_pt[[param_type]], z_zero_ca_pt[[param_type]], z_def_bdry[[param_type]], moments[[param_type]][,'z','mean'])
  pdf(file_name)
  par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
  plot(x = comp_param[[param_type]], y = 100*z_att_pt[[param_type]], lty = 1, 
       bty = 'n', type = 'l',axes = 'F', main = '', xlab = '', ylab = '',col = 'blue', lwd = thickness, cex.lab = cex_plot, 
       cex.axis = cex_plot, ylim = c(y_min,y_max))
  lines(x = comp_param[[param_type]], y = 100*z_zero_drift_pt[[param_type]], col = 'pink', lwd = thickness, lty = 2)
  lines(x = comp_param[[param_type]], y = 100*z_zero_ca_pt[[param_type]], col = 'orange', lwd = thickness, lty = 3)
  lines(x = comp_param[[param_type]], y = 100*moments[[param_type]][,'z','mean'], col = 'green', lwd = thickness, lty = 4)
  lines(x = comp_param[[param_type]], y = 100*z_def_bdry[[param_type]], col = 'red', lwd = thickness, lty = 5)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  axis(2, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(2,text = "debt coverage ratio (%)", line = 2, cex = cex_plot)
  axis(1, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(1,text=param_type, line=2, cex = cex_plot)
  legend("topright", legend = c('attract. point','zero drift','zero current act.','ergodic mean','default point'), 
         col = c('blue', 'pink', 'orange','green','red'), bty="n", lty = 1:5, cex = cex_plot, lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # plot other sensitivities
  mom_names <- c('x','z','spread','cons','ca')
  for (mom in mom_names){
    file_name <- paste0(mom,"_sensi_",param_type,".pdf")
    pdf(file_name)
    par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
    plot(x = comp_param[[param_type]], y = moments[[param_type]][,mom,'mean'], 
         bty = 'n', type = 'l',axes = 'F', main = '', xlab = '', ylab = '',col = 'blue', lwd = thickness, cex.lab = cex_plot, 
         cex.axis = cex_plot)
    if (mom == 'x') lines(x = comp_param[[param_type]], y = att_pt[[param_type]], lwd = thickness, lty = 3, col = 'green')
    if (mom == 'z') lines(x = comp_param[[param_type]], y = z_att_pt[[param_type]], lwd = thickness, lty = 3, col = 'green')
    grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
    axis(2, col="black", lwd = 1, cex.axis = cex_plot)
    mtext(2,text = paste0('average ', mom), line = 2, cex = cex_plot)
    par(new=T)
    plot(x = comp_param[[param_type]], y = moments[[param_type]][,mom,'stdev'], axes=F, xlab="", ylab="", 
         type="l", main="", col='black', lty=2, lwd = thickness)
    axis(1, col="black", lwd = 1, cex.axis = cex_plot)
    mtext(1,text=param_type, line=2, cex = cex_plot)
    axis(4, col="black", lwd = 1, cex.axis = cex_plot)
    mtext(4,text = paste0('stdev ', mom), line = 2, cex = cex_plot)
    if (mom == 'x'){
      legend("bottomright", legend = c('average', 'attraction pt','stdev'), 
             col = c('blue', 'green','black'), bty="n", lty = c(1,3,2), cex = cex_plot, lwd = thickness) 
    } else if (mom == 'z'){
      legend("bottomright", legend = c('average', 'attraction pt','stdev'), 
             col = c('blue', 'green','black'), bty="n", lty = c(1,3,2), cex = cex_plot, lwd = thickness) 
    } else {
      legend("bottomright", legend = c('average', 'stdev'), 
             col = c('blue', 'black'), bty="n", lty = c(1,2), cex = cex_plot, lwd = thickness)
    }
    dev.off()
  }
}

#=============================
# Markov Switching Equilibrium
#=============================

# create output directory for pictures for comparative statics
setwd(output_dir)
if (!file.exists('markov_switching')) dir.create(file.path(output_dir, 'markov_switching'))

# store output in appropriate folder
analysis_dir    <- file.path(output_dir, 'markov_switching')
setwd(analysis_dir)

# parameters
param_4              <- param_0
param_4['theta']     <- .000001
param_4['alpha']     <- .000001
param_4['nu']        <- 0
param_4['mu']        <- 0
param_4['sigma']     <- .2
param_4['r']         <- .05
param_4['delta']     <- .1
param_4['hat_delta'] <- param_4['r']
param_4['kappa']     <- param_4['r']
param_4['m']         <- 1/10

# introduce switching intensities
param_4['lambda_c'] <- .2      # constrained to unconstrained switching intensity
param_4['lambda_u'] <- .2      # unconstrained to constrained switching intensity

# solve for equilibrium cutoff
x_cutoff <- no_commit_cutoff_fn(param_4)

# construct grid
x_min    <- 0.0001
x_max    <- x_cutoff
nb_pts   <- 1000

# calculate all required equilibrium functions
x      <- seq(from = x_min, to = x_max, length.out = nb_pts)
d_u    <- markov_switching_debt_px_u_fn(param_4, x)
d_c    <- markov_switching_debt_px_c_fn(param_4, x)
g_c    <- markov_switching_g_fn(param_4, x)
g_u    <- no_commit_g_fn(param_4, x)

# generate plots in proper directory
setwd(analysis_dir)

# plot debt price
pdf("debt_price_markov_switching.pdf")
plot(x = x, ylim = c(0,1), lwd = thickness, y = d_u, type='l', col = 'blue',
     bty = 'n', xlab = "debt-to-income x", ylab = "debt price", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = x, y = d_c, lwd = thickness, col = 'red', lty = 2)
legend("bottomleft",legend=c(expression("d"[u]*"(x) = d(x)"), 
                             expression("d"[c]*"(x)")), lwd = thickness, 
       col=c('blue','red'), lty= 1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot debt price in eps format
setEPS()
postscript("debt_price_markov_switching.eps")
plot(x = x, ylim = c(0,1), lwd = thickness, y = d_u, type='l', col = 'blue',
     bty = 'n', xlab = "debt-to-income x", ylab = "debt price", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = x, y = d_c, lwd = thickness, col = 'red', lty = 2)
legend("bottomleft",legend=c(expression("d"[u]*"(x) = d(x)"), 
                             expression("d"[c]*"(x)")), lwd = thickness, 
       col=c('blue','red'), lty= 1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot issuance function
pdf("g_markov_switching.pdf")
plot(x = x, lwd = thickness, y = g_u, type='l', cex.lab = cex_plot, lty = 1, 
     cex.axis = cex_plot, col = 'blue', bty = 'n', xlab = "debt-to-income x", 
     ylab = "issuance policy", ylim = c(min(c(g_c,g_u),na.rm = T),10))
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = x, y = g_c, lwd = thickness, col = 'red', lty = 2)
legend("bottomleft",
       legend=c(expression(g*"(x)"), expression(g[u]*"(x)")), lwd = thickness, 
       col=c('blue','red'), lty= 1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot issuance function in eps format
setEPS()
postscript("g_markov_switching.eps")
plot(x = x, lwd = thickness, y = g_u, type='l', cex.lab = cex_plot, lty = 1, 
     cex.axis = cex_plot, col = 'blue', bty = 'n', xlab = "debt-to-income x", 
     ylab = "issuance policy", ylim = c(min(c(g_c,g_u),na.rm = T),10))
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = x, y = g_c, lwd = thickness, col = 'red', lty = 2)
legend("bottomleft",
       legend=c(expression(g*"(x)"), expression(g[u]*"(x)")), lwd = thickness, 
       col=c('blue','red'), lty= 1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# generate csv file
setwd(analysis_dir)
output <- data.frame('d_c' = d_c, 'd_u' = d_u, 'g_c' = g_c, 'g_u' = g_u)
write.table(output, file = "markov_switching_base_case.csv", eol = "\n", na = "NA", fileEncoding = "", row.names = F)

#==============================================
# Debt Ceiling Policies
# Constraint above bar_x* (smooth eqm)
# "Loose Debt Ceiling"
#==============================================

# create output directory for pictures for comparative statics
setwd(output_dir)
if (!file.exists('debt_ceiling')) dir.create(file.path(output_dir, 'debt_ceiling'))

# store output in appropriate folder
analysis_dir    <- file.path(output_dir, 'debt_ceiling')
setwd(analysis_dir)

# parameters
param_5              <- param_4

# solve for equilibrium cutoff
x_cutoff     <- no_commit_cutoff_fn(param_5)
x_star_min   <- issuance_constraint_ratio_limit_fn(param_5)*x_cutoff
x_low_star   <- intermed_x_star_constraint_fn(param_5)

# check whether the intermediate region exists
x_low_star < x_star_min - .01

# set an arbitrary x* > x_star_min, so that smooth equilibrium exists
x_star       <- x_cutoff*(.05)+x_star_min*.95

# construct grid
x_min    <- 0.0001
x_max    <- x_cutoff
nb_pts   <- 1000

# calculate all required equilibrium functions
x       <- seq(from = x_min, to = x_max, length.out = nb_pts)
d_c     <- constrained_smooth_debt_px_fn(param_5, x_star, x)
d_u     <- no_commit_debt_px_fn(param_5, x)
hat_v_c <- constrained_hat_value_fn(param_5, x_star, num_param_0)
hat_v_c0<- hat_v_c[['hat_v0']]
hat_v_c <- hat_v_c[['hat_v']](x)
v_u     <- no_commit_value_fn(param_5, x)

# plot debt price
setwd(analysis_dir)
pdf("debt_price_smooth_constrained.pdf")
plot(x = x, y = d_u, ylim = c(0,1), xlim = c(0.8*x_star, x_max), type='l', col = 'blue', lty = 1, 
     lwd = thickness, bty = 'n', xlab = "debt-to-income x", ylab = "debt price", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = x, y = d_c, lwd = thickness, col = 'red', lty = 2)
legend("topright",legend=c(expression("d(x) (no commit.)"), 
                           expression("d"[c]*"(x) (constr. at x*)")), 
       col=c('blue','red'), lty=1:2, bty = "n", lwd = thickness, cex = cex_plot)
abline(v = x_star, col = 'black', lwd = thickness, lty = 3)
text(x = x_star - .05*(u[2]-u[1]), y = u[4] - .04*(u[4]-u[3]), "x*", cex = cex_plot, col = 'black') 
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot debt price in .eps format
setwd(analysis_dir)
setEPS()
postscript("debt_price_smooth_constrained.eps")
plot(x = x, y = d_u, ylim = c(0,1), xlim = c(0.8*x_star, x_max), type='l', col = 'blue', lty = 1, 
     lwd = thickness, bty = 'n', xlab = "debt-to-income x", ylab = "debt price", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = x, y = d_c, lwd = thickness, col = 'red', lty = 2)
legend("topright",legend=c(expression("d(x) (no commit.)"), 
                           expression("d"[c]*"(x) (constr. at x*)")), 
       col=c('blue','red'), lty=1:2, bty = "n", lwd = thickness, cex = cex_plot)
abline(v = x_star, col = 'black', lwd = thickness, lty = 3)
text(x = x_star - .05*(u[2]-u[1]), y = u[4] - .04*(u[4]-u[3]), "x*", cex = cex_plot, col = 'black') 
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# export
setwd(analysis_dir)
output <- data.frame('x' = x, 'd_c' = d_c, 'd_u' = d_u)
write.table(output, file = "debt_ceiling_smooth_MPE.csv", eol = "\n", na = "NA", fileEncoding = "", row.names = F)

#===========================================================
# Debt Ceiling Policies
# Constraint below bar_x* (reflecting eqm with jump @ x = 0)
# "Tight Debt Ceiling"
#===========================================================

# set x_star to be below x_low_star so that equilibrium is for sure a jump eq.
x_star      <- x_low_star*.75

# solve for equilibrium cutoff
x_cutoff    <- constrained_reflecting_cutoff_fn(param_5,x_star)

# construct grid
x_min    <- .0001
x_max    <- x_cutoff
nb_pts   <- 1000

# calculate all required equilibrium functions
x      <- seq(from = x_min, to = x_max, length.out = nb_pts)
d_c    <- constrained_reflecting_debt_px_fn(param_5, x_star, x)
v_c    <- constrained_reflecting_value_fn(param_5, x_star, x)
d_u    <- no_commit_debt_px_fn(param_5, x)
v_u    <- no_commit_value_fn(param_5, x)

# verifications
eps <- .000001

# d'(x*) = 0
(constrained_reflecting_debt_px_fn(param_5, x_star, x_star+eps)-
    constrained_reflecting_debt_px_fn(param_5, x_star, x_star-eps))/(2*eps)

# v'(x_bar) = 0
(constrained_reflecting_value_fn(param_5, x_star, x_cutoff)-
    constrained_reflecting_value_fn(param_5, x_star, x_cutoff-eps))/eps

# d(x_bar) = 0
constrained_reflecting_debt_px_fn(param_5, x_star, x_cutoff-eps)

# v'(x*) + d(x*) = 0
(constrained_reflecting_value_fn(param_5, x_star, x_star+eps)-
    constrained_reflecting_value_fn(param_5, x_star, x_star-eps))/(2*eps) + 
  constrained_reflecting_debt_px_fn(param_5, x_star, x_star)

# plot debt price
pdf("debt_price_reflecting_constrained.pdf")
plot(x = x[x <= no_commit_cutoff_fn(param_5)], y = d_u[x <= no_commit_cutoff_fn(param_5)], 
     xlim = c(0,x_max), ylim = c(0,1), lwd = thickness, type='l', col = 'blue', lty = 1,
     bty = 'n', xlab = "debt-to-income x", ylab = "debt price", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = x, y = d_c, lwd = thickness, col = 'red', lty = 2)
legend("bottomleft",legend=c(expression("d(x) (no commit.)"), 
                             expression("d"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
abline(v = x_star, col = 'black', lwd = thickness, lty = 3)
text(x = x_star - .05*(u[2]-u[1]), y = u[3] + .5*(u[4]-u[3]), 
     "x*", cex = cex_plot, col = 'black') 
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot debt price in .eps format
setEPS()
postscript("debt_price_reflecting_constrained.eps")
plot(x = x[x <= no_commit_cutoff_fn(param_5)], y = d_u[x <= no_commit_cutoff_fn(param_5)], 
     xlim = c(0,x_max), ylim = c(0,1), lwd = thickness, type='l', col = 'blue', lty = 1,
     bty = 'n', xlab = "debt-to-income x", ylab = "debt price", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = x, y = d_c, lwd = thickness, col = 'red', lty = 2)
legend("bottomleft",legend=c(expression("d(x) (no commit.)"), 
                             expression("d"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
abline(v = x_star, col = 'black', lwd = thickness, lty = 3)
text(x = x_star - .05*(u[2]-u[1]), y = u[3] + .5*(u[4]-u[3]), 
     "x*", cex = cex_plot, col = 'black') 
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot value function
pdf("value_function_reflecting_constrained.pdf")
plot(x = x[x <= no_commit_cutoff_fn(param_5)], y = v_u[x <= no_commit_cutoff_fn(param_5)], 
     ylim = c(min(v_u,v_c),max(v_u,v_c)), xlim = c(0,x_max),  
     type='l', col = 'blue', lwd = thickness, lty = 1, 
     bty = 'n', xlab = "debt-to-income x", ylab = "value function", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = x, y = v_c, lwd = thickness, col = 'red', lty = 2)
legend("topright",legend=c(expression("v(x) (no commit.)"), 
                           expression("v"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
abline(v=x_star, lty = 3, lwd = thickness, col ='black')
text(x = x_star + .05*(u[2]-u[1]), y = u[3] + .04*(u[4]-u[3]), "x*", cex = cex_plot)
arrows(x0 = u[1] + .02*(u[2]-u[1]), y0 = no_commit_value_fn(param_5, 0.000001),
       x1 = u[1] + .02*(u[2]-u[1]), 
       y1 = constrained_reflecting_value_fn(param_5, x_star, 0.000001), 
       angle = 10, code = 3, col = 'black', lty = 1, lwd = thickness, length = .1)
text(x = u[1] + .02*(u[2]-u[1]), y = .5*no_commit_value_fn(param_5, 0.000001) + 
       .5*constrained_reflecting_value_fn(param_5, x_star, 0.000001), "welfare gain", cex = cex_plot, pos = 4)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot value function in .eps format
setEPS()
postscript("value_function_reflecting_constrained.eps")
plot(x = x[x <= no_commit_cutoff_fn(param_5)], y = v_u[x <= no_commit_cutoff_fn(param_5)], 
     ylim = c(min(v_u,v_c),max(v_u,v_c)), xlim = c(0,x_max),  
     type='l', col = 'blue', lwd = thickness, lty = 1, 
     bty = 'n', xlab = "debt-to-income x", ylab = "value function", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = x, y = v_c, lwd = thickness, col = 'red', lty = 2)
legend("topright",legend=c(expression("v(x) (no commit.)"), 
                           expression("v"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
abline(v=x_star, lty = 3, lwd = thickness, col ='black')
text(x = x_star + .05*(u[2]-u[1]), y = u[3] + .04*(u[4]-u[3]), "x*", cex = cex_plot)
arrows(x0 = u[1] + .02*(u[2]-u[1]), y0 = no_commit_value_fn(param_5, 0.000001),
       x1 = u[1] + .02*(u[2]-u[1]), 
       y1 = constrained_reflecting_value_fn(param_5, x_star, 0.000001), 
       angle = 10, code = 3, col = 'black', lty = 1, lwd = thickness, length = .1)
text(x = u[1] + .02*(u[2]-u[1]), y = .5*no_commit_value_fn(param_5, 0.000001) + 
       .5*constrained_reflecting_value_fn(param_5, x_star, 0.000001), "welfare gain", cex = cex_plot, pos = 4)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# export
setwd(analysis_dir)
output <- data.frame('x' = x, 'd_c' = d_c, 'd_u' = d_u, 'v_c' = v_c, 'v_u' = v_u)
write.table(output, file = "debt_ceiling_reflecting_MPE.csv", eol = "\n", na = "NA", fileEncoding = "", row.names = F)

#========================================================
# Debt Ceiling Policies
# Constraint below x* (reflecting eqm with jump @ x = 0)
# Range of x* below limiting value
#========================================================

# set appropriate folder
setwd(analysis_dir)

# vector of x* for which generate plots
x_star    <- seq(from = x_low_star*.95, to = x_star_min*.999, length.out = 20)
x_cutoff  <- array(NA,length(x_star))
x_hat     <- array(NA,length(x_star))
d_x_star  <- array(NA,length(x_star))

# create relevant output directory
setwd(analysis_dir)
if (!file.exists("v_reflecting")) dir.create(file.path(analysis_dir, "v_reflecting"))
if (!file.exists("v_differential_reflecting")) dir.create(file.path(analysis_dir, "v_differential_reflecting"))
if (!file.exists("d_reflecting")) dir.create(file.path(analysis_dir, "d_reflecting"))

for (i in (1:length(x_star))){
  x_st <- x_star[i]
  
  # solve for equilibrium cutoff
  x_cutoff[i]  <- constrained_reflecting_cutoff_fn(param_5,x_st)
  x_hat[i]     <- constrained_reflecting_jump_cutoff_fn2(param_5,x_st)
  
  # construct grid
  x_min    <- .0001
  x_max    <- x_cutoff[i]
  nb_pts   <- 1000
  
  # calculate all required equilibrium functions
  x      <- seq(from = x_min, to = x_max, length.out = nb_pts)
  d_c    <- constrained_reflecting_debt_px_fn(param_5, x_st, x)
  v_c    <- constrained_reflecting_value_fn(param_5, x_st, x)
  d_u    <- no_commit_debt_px_fn(param_5, x)
  v_u    <- no_commit_value_fn(param_5, x)
  
  # compute debt price at x*, assuming reflection
  d_x_star[i]  <- constrained_reflecting_debt_px_fn(param_5, x_st, x_st)
  
  # plot debt price
  setwd(file.path(analysis_dir, "d_reflecting"))
  pdf(paste("debt_price_reflecting_constrained_",i,".pdf",sep=""))
  plot(x = x, y = d_u, ylim = c(0,1), xlim = c(0,1.2*no_commit_cutoff_fn(param_5)), 
       lwd = thickness, type='l', col = 'blue',
       bty = 'n', xlab = "debt-to-income x", ylab = "debt price", cex.lab = cex_plot, cex.axis = cex_plot)
  lines(x = x, y = d_c, lty = 2, lwd = thickness, col = 'red')
  legend("bottomleft",legend=c(expression("d(x) (no commit.)"), 
                               expression("d"[c]*"(x) (constr.)")), lwd = thickness, 
         col=c('blue','red'), lty=1, bty = "n",cex = cex_plot)
  abline(v=no_commit_cutoff_fn(param_5), lty = 3, lwd = thickness, col ='blue')
  abline(v=x_cutoff[i], lty = 3, lwd = thickness, col ='red')
  abline(v=x_st, lty = 3, lwd = thickness, col ='black')
  dev.off()
  
  # plot value function
  setwd(file.path(analysis_dir, "v_reflecting"))
  pdf(paste("value_function_reflecting_constrained_",i,".pdf",sep=""))
  plot(x = x[x <= no_commit_cutoff_fn(param_5)], y = v_u[x <= no_commit_cutoff_fn(param_5)], 
       ylim = c(min(v_u,v_c),max(v_u,v_c)), 
       type='l', col = 'blue', lwd = thickness, xlim = c(0,1.2*no_commit_cutoff_fn(param_5)), 
       bty = 'n', xlab = "debt-to-income x", ylab = "value function", cex.lab = cex_plot, cex.axis = cex_plot)
  u <- par('usr')
  lines(x = x, y = v_c, lty = 2, lwd = thickness, col = 'red')
  legend("topright",legend=c(expression("v(x) (no commit.)"), 
                             expression("v"[c]*"(x) (constr.)")), lwd = thickness, 
         col=c('blue','red'), lty=1, bty = "n",cex = cex_plot)
  abline(v=no_commit_cutoff_fn(param_5), lty = 3, lwd = thickness, col ='blue')
  abline(v=x_cutoff[i], lty = 3, lwd = thickness, col ='red')
  abline(v=x_st, lty = 3, lwd = thickness, col ='black')
  text(x = x_st + .05*(u[2]-u[1]), y = u[3] + .04*(u[4]-u[3]), "x*", cex = cex_plot)
  dev.off()
  
  # plot value function differential
  setwd(file.path(analysis_dir, "v_differential_reflecting"))
  pdf(paste("value_differential_reflecting_constrained_",i,".pdf",sep=""))
  plot(x = x[x <= no_commit_cutoff_fn(param_5)], y = (v_c-v_u)[x <= no_commit_cutoff_fn(param_5)], 
       type='l', col = 'blue', lwd = thickness, xlim = c(0,1.2*no_commit_cutoff_fn(param_5)), 
       bty = 'n', xlab = "debt-to-income x", ylab = expression("v"[c]*"(x) - v(x)"), cex.lab = cex_plot, cex.axis = cex_plot)
  u <- par('usr')
  abline(v=no_commit_cutoff_fn(param_5), lty = 2, lwd = thickness, col ='blue')
  abline(v=x_cutoff[i], lty = 2, lwd = thickness, col ='red')
  abline(v=x_st, lty = 2, lwd = thickness, col ='black')
  if (x_hat[i] > 0) {
    abline(v=x_hat[i], lty = 2, lwd = thickness, col ='lightblue')
    text(x = x_hat[i] - .05*(u[2]-u[1]), y = u[3] + .04*(u[4]-u[3]), expression(hat(x)), cex = cex_plot)
  }
  text(x = x_st + .05*(u[2]-u[1]), y = u[3] + .04*(u[4]-u[3]), "x*", cex = cex_plot)
  dev.off()
}

# plot default boundary vs. x*
setwd(analysis_dir)
pdf("default_boundary_vs_x_star.pdf")
plot(x = x_star, y = x_cutoff, type='l', col = 'blue', lwd = thickness, 
     ylim = c(no_commit_cutoff_fn(param_5),max(x_cutoff)),
     bty = 'n', xlab = "x*", ylab = expression(bar(x)), cex.lab = cex_plot, cex.axis = cex_plot)
abline(h = no_commit_cutoff_fn(param_5), lty = 2, lwd = thickness, col = 'red')
dev.off()

# plot to verify that there is never a smooth region between 0 and x* -- in plot, need to make sure that x_hat is not in [0,x*]
pdf("x_hat_vs_x_star.pdf")
plot(x = x_star, y = x_hat, type='l', col = 'blue', lwd = thickness, ylim = c(min(x_star,x_hat), max(x_star,x_hat)),
     bty = 'n', xlab = "x*", ylab = '', cex.lab = cex_plot, cex.axis = cex_plot)
u <- par('usr')
title(ylab = expression(hat(x)), line = 2, cex.lab = cex_plot)
lines(x = x_star, y = x_star, lty = 2, lwd = thickness, col = 'red')
abline(h = 0, lty = 2, lwd = thickness, col = 'purple')
abline(v = x_star_min, lty = 2, lwd = thickness, col = 'green')
text(x = x_star_min - .05*(u[2]-u[1]), y = u[3] + .04*(u[4]-u[3]), expression(bar(x)~"*"), cex = cex_plot)
dev.off()

# plot debt price at x* -- need to make sure debt price at x* > (kappa+m)/(delta+m), otherwise there is a smooth region near x = 0
d_risk_free_delta <- (param_5['kappa']+param_5['m'])/(param_5['delta']+param_5['m'])
pdf("d_at_x_star.pdf")
plot(x = x_star, y = d_x_star, type='l', col = 'blue', lwd = thickness, ylim = c(min(d_risk_free_delta,d_x_star),max(d_risk_free_delta,d_x_star)),
     bty = 'n', xlab = "x*", ylab = '', cex.lab = cex_plot, cex.axis = cex_plot)
title(ylab = expression(d("x*")), line = 2, cex.lab = cex_plot)
abline(h = d_risk_free_delta, lty = 2, lwd = thickness, col = 'purple')
dev.off()

#=================================================
# Debt Ceiling Policies
# Constraint above bar_x* (smooth eqm)
# Government vs. citizen's welfare for range of x*
# Computations for all values of x
#=================================================

# set appropriate folder
setwd(analysis_dir)

# solve for equilibrium cutoff
x_cutoff     <- no_commit_cutoff_fn(param_5)
x_star_min   <- issuance_constraint_ratio_limit_fn(param_5)*x_cutoff
x_low_star   <- intermed_x_star_constraint_fn(param_5)

# take 2 different hat_delta values
hat_delta_1  <- .5*param_5['r'] + .5*param_5['delta']
hat_delta_2  <- param_5['r']

# construct vector of x_star -- different number of points for different regions
nb_x_star1   <- 50  # number of different policies to compute for x* < x_low_star
nb_x_star2   <- 0   # number of different policies to compute for x_star_min > x* > x_low_star
nb_x_star3   <- 50  # number of different policies to compute for x_star_min > x* > x_low_star
x_eval       <- seq(from = .01, to = 5*x_cutoff/6, length.out = 6) # number of different debt-to-income at which estimating welfare gains/losses
nb_x         <- length(x_eval) 
x_star       <- c(seq(from = x_low_star*.01, to = x_low_star*.99, length.out = nb_x_star1), 
                  # seq(from = x_star_min*(.025)+x_low_star*.975, to = x_low_star*(.025)+x_star_min*.975, length.out = nb_x_star2),
                  seq(from = x_cutoff*(.01)+x_star_min*.99, to = x_cutoff*(.99)+x_star_min*.01, length.out = nb_x_star3))
nb_x_star    <- length(x_star)

# store value function outputs
hat_v1           <- array(NA,dim = c(nb_x_star,nb_x*2))
colnames(hat_v1) <- c(paste0('welfare_u_',(1:nb_x)),
                      paste0('welfare_c_',(1:nb_x)))
hat_v2           <- hat_v1
gov_v            <- hat_v1

# store welfare losses
welfare_chg_1           <- array(NA,dim = c(nb_x_star,nb_x))
colnames(welfare_chg_1) <- paste0('welfare_pct_delta_',(1:nb_x))
welfare_chg_2           <- welfare_chg_1
welfare_chg_gov         <- welfare_chg_1

# store boundaries
bdries                  <- array(NA,dim = c(nb_x_star,2))
colnames(bdries)        <- c('x_bar','x_hat')

# store hat_value functions for checking purposes
x_list      <- list()
hat_v1_list <- list()
hat_v2_list <- list()

# create directories that are needed
for (i in (1:nb_x_star)){
  # print run number
  print(paste("run ",i,sep=""))
  
  # set constraint
  x_st <- x_star[i]
  
  # store x_list
  if (x_star_min>x_st){
    bdries[i,'x_bar'] <- constrained_reflecting_cutoff_fn(param_5,x_st)
    if (x_st>x_low_star){
      # in case intermediate region exists, i.e. "moderate debt ceiling"
      bdries[i,'x_hat'] <- constrained_reflecting_jump_cutoff_fn2(param_5,x_st)
    }
  } else {
    bdries[i,'x_bar'] <- no_commit_cutoff_fn(param_5)
  }
  
  # store x vector
  x_list[[i]] <- seq(from = 0.01, to = bdries[i,'x_bar'], length.out = 1000)
  
  # fill the table for government value function
  for (j in (1:nb_x)){
    gov_v[i,paste0('welfare_u_',j)] <- no_commit_value_fn(param_5,x_eval[j])
    if (x_star_min>x_st){
      # case where equilibrium is a reflecting equilibrium
      gov_v[i,paste0('welfare_c_',j)] <- constrained_reflecting_value_fn(param_5, x_st, x_eval[j])
    } else {
      # case where equilibrium is a smooth equilibrium
      gov_v[i,paste0('welfare_c_',j)] <- no_commit_value_fn(param_5,x_eval[j])
    }
    # compute welfare gain/loss vs. unconstrained
    welfare_chg_gov[i,paste0('welfare_pct_delta_',j)] <- gov_v[i,paste0('welfare_c_',j)]/
      gov_v[i,paste0('welfare_u_',j)] - 1
  }
  
  # calculate hat_v_u and hat_v_c for hat_delta_1
  param_5['hat_delta']  <- hat_delta_1
  
  # compute hat_values_1 for (a) unconstrained eq, (b) constrained eq
  hat_v_u               <- no_commit_hat_value_fn(param_5, num_param_0)[['hat_v']]
  hat_v_c               <- constrained_hat_value_fn(param = param_5, x_star = x_st, num_param = num_param_0)[['hat_v']]
  hat_v1_list[[i]]      <- hat_v_c(x_list[[i]])
  
  # fill the table
  for (j in (1:nb_x)){
    hat_v1[i,paste0('welfare_u_',j)] <- hat_v_u(x_eval[j])
    hat_v1[i,paste0('welfare_c_',j)] <- hat_v_c(x_eval[j])
    # compute welfare loss at 2 different levels of indebtedness
    welfare_chg_1[i,paste0('welfare_pct_delta_',j)] <- hat_v1[i,paste0('welfare_c_',j)]/
      hat_v1[i,paste0('welfare_u_',j)]- 1
  }
  
  # calculate hat_v_u and hat_v_c for hat_delta_2
  param_5['hat_delta']  <- hat_delta_2
  
  # compute hat_values_2 for (a) unconstrained eq, (b) constrained eq
  hat_v_u               <- no_commit_hat_value_fn(param_5, num_param_0)[['hat_v']]
  hat_v_c               <- constrained_hat_value_fn(param_5, x_st, num_param_0)[['hat_v']]
  hat_v2_list[[i]]      <- hat_v_c(x_list[[i]])
  
  # fill the table
  for (j in (1:nb_x)){
    hat_v2[i,paste0('welfare_u_',j)] <- hat_v_u(x_eval[j])
    hat_v2[i,paste0('welfare_c_',j)] <- hat_v_c(x_eval[j])
    # compute welfare loss at 2 different levels of indebtedness
    welfare_chg_2[i,paste0('welfare_pct_delta_',j)] <- hat_v2[i,paste0('welfare_c_',j)]/
      hat_v2[i,paste0('welfare_u_',j)]- 1
  }
}

# plot welfare changes (x_star in horizontal axis)
setwd(analysis_dir)
if (!file.exists("welfare_analysis")) dir.create(file.path(analysis_dir, "welfare_analysis"))
setwd(file.path(analysis_dir, "welfare_analysis"))
for (j in (1:nb_x)){
  y_min <- 100*min(welfare_chg_1[,j], welfare_chg_2[,j],welfare_chg_gov[,j])
  y_max <- 100*max(welfare_chg_1[,j], welfare_chg_2[,j],welfare_chg_gov[,j])
  pdf(paste0("welfare_change_debt_ceiling_",j,".pdf"))
  plot(x = x_star[1:nb_x_star1], y = 100*welfare_chg_1[1:nb_x_star1,j], lty = 2, 
       type='l', xlab = 'x*', ylab = 'welfare gains (%)',
       col="red", lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n',
       ylim = c(y_min,y_max), xlim = c(0,max(x_star)))
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  lines(x = x_star[(nb_x_star1+nb_x_star2+1):nb_x_star], y = 100*welfare_chg_1[(nb_x_star1+nb_x_star2+1):nb_x_star,j], type='l', col="red", lwd = thickness, lty = 2)
  lines(x = x_star[1:nb_x_star1], y = 100*welfare_chg_2[1:nb_x_star1,j], type='l', col="purple", lwd = thickness, lty = 3)
  lines(x = x_star[(nb_x_star1+nb_x_star2+1):nb_x_star], y = 100*welfare_chg_2[(nb_x_star1+nb_x_star2+1):nb_x_star,j], type='l', col="purple", lwd = thickness, lty = 3)
  lines(x = x_star, y = 100*welfare_chg_gov[,j], type='l', col="blue", lwd = thickness, lty = 1)
  u <- par('usr')
  abline(v = x_star_min, lty = 4, lwd = thickness, col ='black')
  text(x = x_star_min - .02*(u[2]-u[1]), y = u[3] + .05*(u[4]-u[3]), expression(bar(x)*'*'), cex = cex_plot, col = 'black') 
  text(x = x_star_min - .5*(x_star_min-u[1]), y = u[3] + .2*(u[4]-u[3]), "reflecting eqm", cex = cex_plot, col = 'black') 
  text(x = x_star_min + .5*(u[2]-x_star_min), y = u[3] + .2*(u[4]-u[3]), "smooth eqm", cex = cex_plot, col = 'black') 
  legend("topright",legend=expression('government', 'citizens ('*delta*'>'*hat(delta)[1]*'>r)', 
                                      'citizens ('*hat(delta)[2]*'=r)'), lwd = thickness, 
         col=c('blue','red','purple'), lty=1:3, bty = "n",cex = cex_plot)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # export
  output <- data.frame('x_star' = x_star, 'welfare_change_1' = welfare_chg_1[,j],
                       'welfare_change_2' = welfare_chg_2[,j], 'welfare_chg_gov' = welfare_chg_gov[,j])
  write.table(output, file = paste0("welfare_at_x_",j,".csv"), eol = "\n", na = "NA", fileEncoding = "", row.names = F)
}

# print eps file for paper
j <- 1
y_min <- 100*min(welfare_chg_1[,j], welfare_chg_2[,j],welfare_chg_gov[,j])
y_max <- 100*max(welfare_chg_1[,j], welfare_chg_2[,j],welfare_chg_gov[,j])
setEPS()
postscript("welfare_change_debt_ceiling_1.eps")
plot(x = x_star[1:nb_x_star1], y = 100*welfare_chg_1[1:nb_x_star1,j], lty = 2, 
     type='l', xlab = 'x*', ylab = 'welfare gains (%)',
     col="red", lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n',
     ylim = c(y_min,y_max), xlim = c(0,max(x_star)))
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = x_star[(nb_x_star1+nb_x_star2+1):nb_x_star], y = 100*welfare_chg_1[(nb_x_star1+nb_x_star2+1):nb_x_star,j], type='l', col="red", lwd = thickness, lty = 2)
lines(x = x_star[1:nb_x_star1], y = 100*welfare_chg_2[1:nb_x_star1,j], type='l', col="purple", lwd = thickness, lty = 3)
lines(x = x_star[(nb_x_star1+nb_x_star2+1):nb_x_star], y = 100*welfare_chg_2[(nb_x_star1+nb_x_star2+1):nb_x_star,j], type='l', col="purple", lwd = thickness, lty = 3)
lines(x = x_star, y = 100*welfare_chg_gov[,j], type='l', col="blue", lwd = thickness, lty = 1)
u <- par('usr')
abline(v = x_star_min, lty = 4, lwd = thickness, col ='black')
text(x = x_star_min - .02*(u[2]-u[1]), y = u[3] + .05*(u[4]-u[3]), expression(bar(x)*'*'), cex = cex_plot, col = 'black') 
text(x = x_star_min - .5*(x_star_min-u[1]), y = u[3] + .2*(u[4]-u[3]), "reflecting eqm", cex = cex_plot, col = 'black') 
text(x = x_star_min + .5*(u[2]-x_star_min), y = u[3] + .2*(u[4]-u[3]), "smooth eqm", cex = cex_plot, col = 'black') 
legend("topright",legend=expression('government', 'citizens ('*delta*'>'*hat(delta)[1]*'>r)', 
                                    'citizens ('*hat(delta)[2]*'=r)'), lwd = thickness, 
       col=c('blue','red','purple'), lty=1:3, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

#========================================================
# Debt Ceiling Policies
# Constraint below x* (reflecting eqm, smooth @ x = 0)
# Computations and plots
#========================================================

# first, make discount rate of government close to the discount rate of creditors
param_6          <- param_5
param_6['delta'] <- param_6['r'] + .02
param_6['m']     <- 1/20

# compute cutoffs
x_star_min  <- issuance_constraint_ratio_limit_fn(param_6)*no_commit_cutoff_fn(param_6)
x_star      <- x_star_min*.98

# solve for equilibrium cutoff
x_cutoff    <- constrained_reflecting_cutoff_fn(param_6,x_star)
x_low       <- constrained_reflecting_jump_cutoff_fn2(param_6,x_star)

# check that equilibrium indeed has 3 regions
constrained_reflecting_debt_px_fn(param_6, x_star, x_star) < (param_6['kappa']+param_6['m'])/(param_6['delta']+param_6['m'])

# construct grid
x_min    <- .0001
x_max    <- x_cutoff
nb_pts   <- 1000

# calculate all required equilibrium functions
x      <- seq(from = x_min, to = x_max, length.out = nb_pts)
d_c    <- constrained_reflecting_debt_px_fn(param_6, x_star, x)
v_c    <- constrained_reflecting_value_fn(param_6, x_star, x)
d_u    <- no_commit_debt_px_fn(param_6, x)
v_u    <- no_commit_value_fn(param_6, x)

# set directory
setwd(analysis_dir)

# plot debt price
pdf("debt_price_reflecting_constrained2.pdf")
plot(x = x, y = d_u, type='l', col = 'blue', bty = 'n', ylim = c(0,1), lwd = thickness, lty = 1,
     xlab = "debt-to-income x", ylab = "debt price", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = x, y = d_c, lwd = thickness, col = 'red', lty = 2)
abline(v=x_star, lty = 3, lwd = thickness, col ='black')
abline(v=x_low, lty = 4, lwd = thickness, col ='pink')
text(x = x_star + .05*(u[2]-u[1]), y = u[3] + .5*(u[4]-u[3]), "x*", cex = cex_plot)
text(x = x_low - .05*(u[2]-u[1]), y = u[3] + .5*(u[4]-u[3]), 
     expression(hat(x)), cex = cex_plot)
legend("bottomleft",legend=c(expression("d(x) (no commit.)"), 
                             expression("d"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()


# plot debt price in eps format
setEPS()
postscript("debt_price_reflecting_constrained2.eps")
plot(x = x, y = d_u, type='l', col = 'blue', bty = 'n', ylim = c(0,1), lwd = thickness, lty = 1,
     xlab = "debt-to-income x", ylab = "debt price", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = x, y = d_c, lwd = thickness, col = 'red', lty = 2)
abline(v=x_star, lty = 3, lwd = thickness, col ='black')
abline(v=x_low, lty = 4, lwd = thickness, col ='pink')
text(x = x_star + .05*(u[2]-u[1]), y = u[3] + .5*(u[4]-u[3]), "x*", cex = cex_plot)
text(x = x_low - .05*(u[2]-u[1]), y = u[3] + .5*(u[4]-u[3]), 
     expression(hat(x)), cex = cex_plot)
legend("bottomleft",legend=c(expression("d(x) (no commit.)"), 
                             expression("d"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()


# plot value function
pdf("value_function_reflecting_constrained2.pdf")
plot(x = x, y = v_u, ylim = c(min(v_u,v_c),max(v_u,v_c)), type='l', col = 'blue', lwd = thickness, lty = 1,
     bty = 'n', xlab = "debt-to-income x", ylab = "value function", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = x, y = v_c, lwd = thickness, col = 'red', lty = 2)
abline(v=x_low, lty = 3, lwd = thickness, col ='pink')
abline(v=x_star, lty = 4, lwd = thickness, col ='black')
text(x = x_star + .05*(u[2]-u[1]), y = u[3] + .04*(u[4]-u[3]), "x*", cex = cex_plot)
text(x = x_low - .05*(u[2]-u[1]), y = u[3] + .04*(u[4]-u[3]), 
     expression(hat(x)), cex = cex_plot)
legend("topright",legend=c(expression("v(x) (no commit.)"), 
                           expression("v"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot value function differential
diff_v <- v_c - v_u
pdf("value_function_diff_reflecting_constrained2.pdf")
plot(x = x, y = diff_v, ylim = c(min(diff_v,0),max(diff_v,0)), type='l', col = 'blue', lwd = thickness, 
     bty = 'n', xlab = "debt-to-income x", ylab = expression("v"[c]*"(x)" - "v(x)"), cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
abline(h = 0, lwd = thickness, lty = 2, col = 'red')
u <- par('usr')
abline(v=x_low, lty = 3, lwd = thickness, col ='pink')
abline(v=x_star, lty = 4, lwd = thickness, col ='black')
text(x = x_star + .05*(u[2]-u[1]), y = u[3] + .04*(u[4]-u[3]), "x*", cex = cex_plot)
text(x = x_low - .05*(u[2]-u[1]), y = u[3] + .04*(u[4]-u[3]), 
     expression(hat(x)), cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# export
output <- data.frame('x' = x, 'v_c' = v_c, 'v_u' = v_u, 'd_c' = d_c, 'd_u' = d_u)
write.table(output, file = "smooth_and_jump_MPE_illustration.csv", eol = "\n", na = "NA", fileEncoding = "", row.names = F)

#=========================================
# Issuance Cap Equilibrium g(x) < bar_g
# Compute and plot for one particular case
#=========================================

# create output directory for pictures for comparative statics
setwd(output_dir)
if (!file.exists('issuance_cap')) dir.create(file.path(output_dir, 'issuance_cap'))

# store output in appropriate folder
analysis_dir    <- file.path(output_dir, 'issuance_cap')
setwd(analysis_dir)

# model parameters
param_7          <- param_4
param_7['mu']    <- .02
param_7['hat_delta'] <- .5*param_7['r']+.5*param_7['delta']

# parameters related to constrained equilibrium
param_7['bar_g']    <- 1     # maximum issuance rate
param_7['reinject'] <- F     # whether particle is reinjected or not

# compute constrained equilibrium
constrained_eq   <- max_g_eq_fn(param_7, num_param_0)
# constrained_eq2  <- max_g_eq_fn2(param_7, num_param_0)

# set up constrained and unconstrained axis
x_c <- seq(from = .01, to = .99*constrained_eq[["x_bar"]], length.out = 250)
x_u <- seq(from = .01, to = .99*no_commit_cutoff_fn(param_7), length.out = 250)

# calculated constrained and unconstrained functions of interest
d_c   <- constrained_eq[['d']](x_c)
d_u   <- no_commit_debt_px_fn(param_7, x_u)
v_c   <- constrained_eq[['v']](x_c)
v_u   <- no_commit_value_fn(param_7, x_u)
g_c   <- constrained_eq[['g']](x_c)
g_u   <- no_commit_g_fn(param_7, x_u)
c_c   <- constrained_eq[['c']](x_c)
c_u   <- no_commit_c_fn(param_7, x_u)
inc_c <- constrained_eq[['incentive']](x_c)

# plot debt price
pdf("debt_price_issuance_cap.pdf")
plot(x = x_c, xlim = c(0,max(x_c,x_u)), lwd = thickness, y = d_c, type='l', col = 'red', ylim = c(0,max(d_c,d_u)),
     bty = 'n', xlab = "debt-to-income x", ylab = "debt price", cex.lab = cex_plot, cex.axis = cex_plot, lty = 2)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 2, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = x_u, y = d_u, col = 'blue', lwd = thickness, lty = 1)
abline(v = constrained_eq[["x_star"]], col = 'black', lty = 3,lwd = thickness)
text(x = constrained_eq[["x_star"]] + .05*(u[2]-u[1]),  y = u[3] + .04*(u[4]-u[3]), 
     expression("x*"), cex = cex_plot, col = 'black') 
legend("bottomleft",legend=c(expression("d(x) (no commit.)"), 
                             expression("d"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot debt price in eps format
setEPS()
postscript("debt_price_issuance_cap.eps")
plot(x = x_c, xlim = c(0,max(x_c,x_u)), lwd = thickness, y = d_c, type='l', col = 'red', ylim = c(0,max(d_c,d_u)),
     bty = 'n', xlab = "debt-to-income x", ylab = "debt price", cex.lab = cex_plot, cex.axis = cex_plot, lty = 2)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 2, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = x_u, y = d_u, col = 'blue', lwd = thickness, lty = 1)
abline(v = constrained_eq[["x_star"]], col = 'black', lty = 3,lwd = thickness)
text(x = constrained_eq[["x_star"]] + .05*(u[2]-u[1]),  y = u[3] + .04*(u[4]-u[3]), 
     expression("x*"), cex = cex_plot, col = 'black') 
legend("bottomleft",legend=c(expression("d(x) (no commit.)"), 
                             expression("d"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot value function
pdf("value_function_issuance_cap.pdf")
plot(x = x_c, xlim = c(0,max(x_c,x_u)), lwd = thickness, y = v_c, type='l', col = 'red', ylim = c(min(v_c,v_u),max(v_c,v_u)), lty = 2,
     bty = 'n', xlab = "debt-to-income x", ylab = "value function", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = x_u, y = v_u, col = 'blue', lwd = thickness, lty = 1)
u <- par('usr')
abline(v = constrained_eq[["x_star"]], col = 'black', lty = 3, lwd = thickness)
# abline(v = constrained_eq[["x_bar"]], col = 'red', lty = 2,lwd = thickness)
# abline(v = no_commit_cutoff_fn(param_7), col = 'blue', lty = 2,lwd = thickness)
text(x = constrained_eq[["x_star"]] - .05*(u[2]-u[1]),  y = u[3] + .04*(u[4]-u[3]), 
     expression("x*"), cex = cex_plot, col = 'black') 
legend("topright",legend=c(expression("v(x) (no commit.)"), 
                           expression("v"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
arrows(x0 = u[1] + .02*(u[2]-u[1]), y0 = 1/(param_7['delta'] - param_7['mu']),
       x1 = u[1] + .02*(u[2]-u[1]), y1 = constrained_eq[['v']](0), 
       angle = 10, code = 3, col = 'black', lty = 1, lwd = thickness)
text(x = u[1] + .02*(u[2]-u[1]), y = .5/(param_7['delta'] - param_7['mu']) + .5*constrained_eq[['v']](0), 
     "welfare gain", cex = cex_plot, pos = 4)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot value function in eps format
setEPS()
postscript("value_function_issuance_cap.eps")
plot(x = x_c, xlim = c(0,max(x_c,x_u)), lwd = thickness, y = v_c, type='l', col = 'red', ylim = c(min(v_c,v_u),max(v_c,v_u)), lty = 2,
     bty = 'n', xlab = "debt-to-income x", ylab = "value function", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = x_u, y = v_u, col = 'blue', lwd = thickness, lty = 1)
u <- par('usr')
abline(v = constrained_eq[["x_star"]], col = 'black', lty = 3, lwd = thickness)
text(x = constrained_eq[["x_star"]] - .05*(u[2]-u[1]),  y = u[3] + .04*(u[4]-u[3]), 
     expression("x*"), cex = cex_plot, col = 'black') 
legend("topright",legend=c(expression("v(x) (no commit.)"), 
                           expression("v"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
arrows(x0 = u[1] + .02*(u[2]-u[1]), y0 = 1/(param_7['delta'] - param_7['mu']),
       x1 = u[1] + .02*(u[2]-u[1]), y1 = constrained_eq[['v']](0), 
       angle = 10, code = 3, col = 'black', lty = 1, lwd = thickness)
text(x = u[1] + .02*(u[2]-u[1]), y = .5/(param_7['delta'] - param_7['mu']) + .5*constrained_eq[['v']](0), 
     "welfare gain", cex = cex_plot, pos = 4)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot issuance policy
pdf("g_issuance_cap.pdf")
plot(x = x_c, xlim = c(0,max(x_c,x_u)),lwd = thickness, y = g_c, type='l', ylim = c(0,3*param_7['bar_g']), col = 'red', lty = 2, 
     bty = 'n', xlab = "debt-to-income x", ylab = "issuance policy g(x)", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = x_u, y = g_u, col = 'blue', lwd = thickness, lty = 1)
abline(v = constrained_eq[["x_star"]], col = 'black', lty = 3, lwd = thickness)
u <- par('usr')
text(x = constrained_eq[["x_star"]] + .03*(u[2]-u[1]), y = u[3] + .03*(u[4]-u[3]), 
     expression("x*"), cex = cex_plot, col = 'black') 
legend("topleft",legend=c(expression("g(x) (no commit.)"), 
                          expression("g"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot issuance policy in eps format
setEPS()
postscript("g_issuance_cap.eps")
plot(x = x_c, xlim = c(0,max(x_c,x_u)),lwd = thickness, y = g_c, type='l', ylim = c(0,3*param_7['bar_g']), col = 'red', lty = 2, 
     bty = 'n', xlab = "debt-to-income x", ylab = "issuance policy g(x)", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = x_u, y = g_u, col = 'blue', lwd = thickness, lty = 1)
abline(v = constrained_eq[["x_star"]], col = 'black', lty = 3, lwd = thickness)
u <- par('usr')
text(x = constrained_eq[["x_star"]] + .03*(u[2]-u[1]), y = u[3] + .03*(u[4]-u[3]), 
     expression("x*"), cex = cex_plot, col = 'black') 
legend("topleft",legend=c(expression("g(x) (no commit.)"), 
                          expression("g"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot consumption policy
pdf("consumption_issuance_cap.pdf")
plot(x = x_c, xlim = c(0,max(x_c,x_u)),lwd = thickness, y = c_c, type='l', ylim = c(-2,5),
     bty = 'n', xlab = "debt-to-income x", ylab = "consumption-to-income", lty = 2,
     cex.lab = cex_plot, cex.axis = cex_plot, col = 'red')
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = x_u, y = c_u, col = 'blue', lwd = thickness, lty = 1)
u <- par('usr')
abline(v = constrained_eq[["x_star"]], col = 'black', lty = 3, lwd = thickness)
# abline(v = constrained_eq[["x_bar"]], col = 'red', lty = 2,lwd = thickness)
# abline(v = no_commit_cutoff_fn(param_7), col = 'blue', lty = 2,lwd = thickness)
text(x = constrained_eq[["x_star"]] + .03*(u[2]-u[1]), y = u[3] + .03*(u[4]-u[3]), 
     expression("x*"), cex = cex_plot, col = 'black') 
legend("topleft",legend=c(expression("c(x) (no commit.)"), 
                          expression("c"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot consumption policy in eps format
setEPS()
postscript("consumption_issuance_cap.eps")
plot(x = x_c, xlim = c(0,max(x_c,x_u)),lwd = thickness, y = c_c, type='l', ylim = c(-2,5),
     bty = 'n', xlab = "debt-to-income x", ylab = "consumption-to-income", lty = 2,
     cex.lab = cex_plot, cex.axis = cex_plot, col = 'red')
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = x_u, y = c_u, col = 'blue', lwd = thickness, lty = 1)
u <- par('usr')
abline(v = constrained_eq[["x_star"]], col = 'black', lty = 3, lwd = thickness)
text(x = constrained_eq[["x_star"]] + .03*(u[2]-u[1]), y = u[3] + .03*(u[4]-u[3]), 
     expression("x*"), cex = cex_plot, col = 'black') 
legend("topleft",legend=c(expression("c(x) (no commit.)"), 
                          expression("c"[c]*"(x) (constr.)")), lwd = thickness, 
       col=c('blue','red'), lty=1:2, bty = "n",cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# export
output <- data.frame('x_c' = x_c, 'x_u' = x_u, 'v_c' = v_c, 'v_u' = v_u, 'd_c' = d_c, 'd_u' = d_u,
                     'g_c' = g_c, 'g_u' = g_u, 'c_c' = c_c, 'c_u' = c_u)
write.table(output, file = "smooth_and_jump_MPE_illustration.csv", eol = "\n", na = "NA", fileEncoding = "", row.names = F)

#=========================================
# Issuance Cap Equilibrium g(x) < bar_g
# Sensitivity to bar_g
#=========================================

# range of bar_g for which to compute equilibrium
bar_g  <- seq(from = .2, to = 3, length.out = 25)
x_bar  <- array(NA,length(bar_g))
x_star <- array(NA,length(bar_g))

# x at which welfare gains are evaluated
x_eval <- c(0, no_commit_cutoff_fn(param_7)/3, 2*no_commit_cutoff_fn(param_7)/3)
nb_x   <- length(x_eval)

# citizens' welfare -- 2 cases
hat_delta_1  <- .5*param_7['r'] + .5*param_7['delta']
hat_delta_2  <- param_7['r']
hat_v_1      <- array(NA,dim = c(length(bar_g), 2*nb_x))
colnames(hat_v_1) <- c(paste0('welfare_u_',(1:nb_x)),
                       paste0('welfare_c_',(1:nb_x)))
hat_v_2      <- hat_v_1
v_0          <- hat_v_1

# welfare pct change -- 2 cases
welfare_chg_0 <- array(NA,dim = c(length(bar_g), nb_x))
colnames(welfare_chg_0) <- c(paste0('welfare_pct_delta_',(1:nb_x)))
welfare_chg_1 <- welfare_chg_0 
welfare_chg_2 <- welfare_chg_0 

for (i in (1:length(bar_g))){
  print(paste('iteration ',i,sep=''))
  param_7['bar_g'] <- bar_g[i]
  constrained_eq   <- max_g_eq_fn(param_7, num_param_0)
  x_star[i]        <- c(constrained_eq[["x_star"]])
  x_bar[i]         <- c(constrained_eq[["x_bar"]])
  
  # store welfare for the government for all possible x_eval
  for (j in (1:nb_x)){
    v_0[i,paste0('welfare_u_',j)] <- no_commit_value_fn(param_7,x_eval[j])
    v_0[i,paste0('welfare_c_',j)] <- constrained_eq[['v']](x_eval[j])
    welfare_chg_0[i,paste0('welfare_pct_delta_',j)] <- v_0[i,paste0('welfare_c_',j)]/
      v_0[i,paste0('welfare_u_',j)] - 1
  }
  
  # compute hat v1 and hat v2
  x    <- seq(from = .001, to = x_bar[i], length.out = 1000)
  g    <- constrained_eq[['g']](x)
  cons <- constrained_eq[['c']](x)
  
  # compute citizen's welfare when hat delta = hat delta 1 
  param_7['hat_delta'] <- hat_delta_1
  hat_v_c              <- commit_hat_value_fn(param_7, user_spec = 1, x, g, cons, num_param_0)
  hat_v_u              <- no_commit_hat_value_fn(param_7, num_param_0)
  
  # store welfare for citizens for all possible x_eval
  for (j in (1:nb_x)){
    hat_v_1[i,paste0('welfare_c_',j)] <- hat_v_c[['hat_v']](x_eval[j])
    hat_v_1[i,paste0('welfare_u_',j)] <- hat_v_u[['hat_v']](x_eval[j])
    welfare_chg_1[i,paste0('welfare_pct_delta_',j)] <- hat_v_1[i,paste0('welfare_c_',j)]/
      hat_v_1[i,paste0('welfare_u_',j)] - 1
  }
  
  # compute citizen's welfare when hat delta = hat delta 2
  param_7['hat_delta'] <- hat_delta_2
  hat_v_c              <- commit_hat_value_fn(param_7, user_spec = 1, x, g, cons, num_param_0)
  hat_v_u              <- no_commit_hat_value_fn(param_7, num_param_0)
  
  # store welfare for citizens for all possible x_eval
  for (j in (1:nb_x)){
    hat_v_2[i,paste0('welfare_c_',j)] <- hat_v_c[['hat_v']](x_eval[j])
    hat_v_2[i,paste0('welfare_u_',j)] <- hat_v_u[['hat_v']](x_eval[j])
    welfare_chg_2[i,paste0('welfare_pct_delta_',j)] <- hat_v_2[i,paste0('welfare_c_',j)]/
      hat_v_2[i,paste0('welfare_u_',j)] - 1
  }
}

# plot barriers
pdf("barriers_vs_bar_g.pdf")
plot(x = bar_g, lwd = thickness,
     y = x_bar, type='l', col = 'blue', ylim = c(min(x_bar,x_star),max(x_bar,x_star)), lty = 1,
     bty = 'n', xlab = expression(bar(g)), ylab = "barriers", cex.lab = cex_plot, cex.axis = cex_plot)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = bar_g, y = x_star, col = 'green', lwd = thickness, lty = 2)
abline(h = no_commit_cutoff_fn(param_7), col = 'red', lty = 3, lwd = thickness)
legend("bottomleft", legend = c(expression(bar(x)[c]), expression('x*')), 
       col = c('blue', 'green'), bty="n", lty = 1:2, cex = cex_plot,lwd = thickness)
text(y = no_commit_cutoff_fn(param_7) - .05*(u[4]-u[3]), x = u[2] - .05*(u[2]-u[1]), 
     expression(bar(x)), cex = cex_plot, col = 'black') 
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot welfare changes (bar g in horizontal axis)
for (j in (1:nb_x)){
  y_min <- 100*min(welfare_chg_1[,j], welfare_chg_2[,j], welfare_chg_0[,j], 0)
  y_max <- 100*max(welfare_chg_1[,j], welfare_chg_2[,j], welfare_chg_0[,j], 0)
  pdf(paste0("welfare_change_issuance_cap_",j,".pdf"))
  plot(x = bar_g, y = 100*welfare_chg_1[,j], 
       type='l', xlab = expression(bar(g)), ylab = 'welfare gain (%)', lty = 2, 
       col="red", lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n',
       ylim = c(y_min,y_max))
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  lines(x = bar_g, y = 100*welfare_chg_2[,j], type='l', col="purple", lwd = thickness, lty = 3)
  lines(x = bar_g, y = 100*welfare_chg_0[,j], type='l', col="blue", lwd = thickness, lty = 1)
  legend("topright", 
         legend = expression('government', 'citizens ('*delta*'>'*hat(delta)[1]*'>r)', 
                             'citizens ('*hat(delta)[2]*'=r)'), 
         col = c('blue', 'red', 'purple'), bty="n", lty = 1:3, cex = cex_plot, lwd = thickness)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
  
  # export
  output <- data.frame('bar_g' = bar_g, 'x_bar' = x_bar, 'welfare_change_1' = welfare_chg_1[,j],
                       'welfare_change_2' = welfare_chg_2[,j], 'welfare_chg_0' = welfare_chg_0[,j])
  write.table(output, file = paste0("welfare_at_x_",j,".csv"), eol = "\n", na = "NA", fileEncoding = "")
}

# plot welfare changes (bar g in horizontal axis) -- eps format for paper
j <- 1
y_min <- 100*min(welfare_chg_1[,j], welfare_chg_2[,j], welfare_chg_0[,j], 0)
y_max <- 100*max(welfare_chg_1[,j], welfare_chg_2[,j], welfare_chg_0[,j], 0)
setEPS()
postscript("welfare_change_issuance_cap_1.eps")
plot(x = bar_g, y = 100*welfare_chg_1[,j], 
     type='l', xlab = expression(bar(g)), ylab = 'welfare gain (%)', lty = 2, 
     col="red", lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n',
     ylim = c(y_min,y_max))
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = bar_g, y = 100*welfare_chg_2[,j], type='l', col="purple", lwd = thickness, lty = 3)
lines(x = bar_g, y = 100*welfare_chg_0[,j], type='l', col="blue", lwd = thickness, lty = 1)
legend("topright", 
       legend = expression('government', 'citizens ('*delta*'>'*hat(delta)[1]*'>r)', 
                           'citizens ('*hat(delta)[2]*'=r)'), 
       col = c('blue', 'red', 'purple'), bty="n", lty = 1:3, cex = cex_plot, lwd = thickness)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

#================================
# Poisson Issuance Model -- setup
#================================

# create output directory for pictures for Poisson model
setwd(output_dir)
if (!file.exists('poisson_model')) dir.create(file.path(output_dir, 'poisson_model'))

# store output in appropriate folder
analysis_dir    <- file.path(output_dir, 'poisson_model')
setwd(analysis_dir)

# model parameters
param_9            <- param_0
param_9['nu']      <- 0
param_9['lambda']  <- 1
param_9['alpha']   <- .001
param_9['theta']   <- .001
param_9['gamma']   <- 0
param_9['mu']      <- .02
param_9['sigma']   <- .2
param_9['r']       <- .05
param_9['delta']   <- .07
param_9['m']       <- .4

# set the vector of maturities for which we will be computing our equilibria
# set also the number of grid points x at which to evaluate function v
m_vec      <- c(1/15,1/10,1/5)
n_m        <- length(m_vec)
n_x        <- 100

# numerical parameters -- for no commitment model
num_param_2               <- num_param_0
num_param_2['tol']        <- 1E-6
num_param_2['dt']         <- 1
num_param_2['max_iter']   <- 5000
num_param_2['print_flag'] <- F

# numerical parameters -- for commitment model
num_param_3               <- num_param_0
num_param_3['tol']        <- 1E-5
num_param_3['dt']         <- 2
num_param_3['max_iter']   <- 100
num_param_3['print_flag'] <- F

#=======================================================================
# Poisson Issuance Model -- No Commitment -- sensitivity to Poisson rate
#=======================================================================

# range of intensities for which to compute equilibrium
lambda_vec <- c(.005,.01,.025,.05,.1,.25,.5,.75,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,7,8,9,10,1000000)
n_lambda   <- length(lambda_vec)
dt_mat     <- array(c(rep(.05,n_lambda), 
                      rep(.05,n_lambda), 
                      rep(.1,n_lambda)), dim = c(n_lambda,n_m))
x_eval     <- list()
out_no_commit <- array(NA,dim = c(n_lambda,7+n_x,n_m))
dimnames(out_no_commit)[[2]] <- c('dt','m','v0','iter','err','x_bar','x_a',paste0('x_',1:n_x))

# loop over all possible maturities
param_9a   <- param_9
for (m_id in (1:n_m)){
  # update maturity
  param_9a['m'] <- m_vec[m_id]
  
  # update x_eval
  x_eval[[m_id]] <- seq(from = 0.025, to = 1.2*no_commit_cutoff_fn(param_9a), length.out = n_x)
  
  setwd(analysis_dir)
  if (!file.exists(paste0('no_commit_maturity_',m_id))) dir.create(file.path(analysis_dir, paste0('no_commit_maturity_',m_id)))
  setwd(file.path(analysis_dir, paste0('no_commit_maturity_',m_id)))
  output_dir_temp <- file.path(analysis_dir, paste0('no_commit_maturity_',m_id))
  
  # loop over all possible time periods dt -- except the last one
  for (i in (1:(n_lambda-1))){
    print(paste0('completion: ',floor(100*i/n_lambda/n_m),'%'))
    param_9a['lambda'] <- lambda_vec[i]
    num_param_2['dt']  <- dt_mat[i,m_id]
    
    # set appropriate directory
    setwd(output_dir_temp)
    if (!file.exists(paste0('poisson_',i))) dir.create(file.path(output_dir_temp, paste0('poisson_',i)))
    setwd(file.path(output_dir_temp, paste0('poisson_',i)))
    
    # compute solution for the no commitment problem
    no_commit_sol           <- Poisson_no_commit_fn(param = param_9a,
                                                    num_param = num_param_2,
                                                    plotting_param = plot_0)
    
    # store solution of the no-commitment problem
    out_no_commit[i,'dt',m_id]   <- 1/dt_mat[i,m_id]
    out_no_commit[i,'m',m_id]    <- m_vec[m_id]
    out_no_commit[i,'v0',m_id]   <- no_commit_sol[['other_moments']]['v_0']
    out_no_commit[i,'iter',m_id] <- no_commit_sol[['algo_perf']]['nb_iter']
    out_no_commit[i,'err',m_id]  <- no_commit_sol[['algo_perf']]['error']
    out_no_commit[i,'x_bar',m_id]<- no_commit_sol[['other_moments']]['x_bar']
    out_no_commit[i,'x_a',m_id]  <- no_commit_sol[['other_moments']]['x_a']
    out_no_commit[i,paste0('x_',1:n_x),m_id] <- no_commit_sol[['eq_functions']][['v']](x_eval[[m_id]])
    
    # store results
    write.table(out_no_commit[i,c('dt','v0','iter','err','x_bar','x_a'),m_id], file = paste0("poisson_model_no_commit_maturity_",m_id,"_lambda_",i,".csv"), eol = "\n", na = "NA", fileEncoding = "", col.names = F)
  }
  
  # fill no commitment case --- case lambda = +infty
  out_no_commit[n_lambda,'dt',m_id]               <- 0
  out_no_commit[n_lambda,'m',m_id]                <- m_vec[m_id]
  out_no_commit[n_lambda,'v0',m_id]               <- no_commit_value_fn(param = param_9a, 0)
  out_no_commit[n_lambda,'x_bar',m_id]            <- no_commit_cutoff_fn(param = param_9a)
  out_no_commit[n_lambda,'x_a',m_id]              <- no_commit_attraction_pt_fn(param = param_9a)
  out_no_commit[n_lambda,paste0('x_',1:n_x),m_id] <- no_commit_value_fn(param = param_9a, x_eval[[m_id]])
}

# export main results
setwd(analysis_dir)
for (m_id in (1:n_m)){
  write.table(out_no_commit[,c('dt','m','v0','iter','err','x_bar','x_a'),m_id], 
              file = paste0('no_commitment_results_m_id_',m_id,'.csv'), eol = "\n", na = "NA", fileEncoding = "", sep = "|", row.names = FALSE)
}

#=========================================================================
# Poisson Issuance Model -- Full Commitment -- sensitivity to Poisson rate
#=========================================================================

setwd(analysis_dir)
lambda_vec2 <- c(seq(from = .01, to = .1, length.out = 4),
                 seq(from = .2, to = 1, length.out = 9),
                 seq(from = 2, to = 10, length.out = 9),
                 seq(from = 15, to = 50, length.out = 8),
                 seq(from = 75, to = 200, length.out = 6),1000000)
n_lambda2   <- length(lambda_vec2)
x_eval2     <- list()
out_commit  <- array(NA,dim = c(n_lambda2,6+n_x,n_m))
dimnames(out_commit)[[2]] <- c('dt','m','v0','x_bar','x_bar_bar','x_star',paste0('x_',1:n_x))

# set parameters
param_9b   <- param_9
for (m_id in (1:n_m)){
  # update maturity
  param_9b['m']   <- m_vec[m_id]
  
  # update x_eval
  x_eval2[[m_id]] <- seq(from = 0.025, to = 1.2*no_commit_cutoff_fn(param_9b), length.out = n_x)
  
  # set appropriate directory
  setwd(analysis_dir)
  if (!file.exists(paste0('commit_maturity_',m_id))) dir.create(file.path(analysis_dir, paste0('commit_maturity_',m_id)))
  setwd(file.path(analysis_dir, paste0('commit_maturity_',m_id)))
  output_dir_temp <- file.path(analysis_dir, paste0('commit_maturity_',m_id))
  setwd(output_dir_temp)
  
  for (i in (1:(n_lambda2-1))){
    print(paste0('completion: ',floor(100*i/n_lambda2/n_m),'%'))
    param_9b['lambda'] <- lambda_vec2[i]
    
    # set appropriate directory
    setwd(output_dir_temp)
    if (!file.exists(paste0('poisson_',i))) dir.create(file.path(output_dir_temp, paste0('poisson_',i)))
    setwd(file.path(output_dir_temp, paste0('poisson_',i)))
    
    # compute solution for the capital structure commitment problem
    out_commit[i,'x_star',m_id] <- Poisson_commit_optim_leverage_fn2(param = param_9b)
    temp                        <- Poisson_commit_cutoff_fn2(param = param_9b, x_star = out_commit[i,'x_star',m_id])
    
    # plots
    plot_Poisson_commit_fn2(param = param_9b, x_star = out_commit[i,'x_star',m_id], plotting_param = plot_0)
    
    # store results
    out_commit[i,'m',m_id]                <- m_vec[m_id]
    out_commit[i,'x_bar',m_id]            <- temp['x_bar']
    out_commit[i,'x_bar_bar',m_id]        <- temp['x_bar_bar']
    out_commit[i,'v0',m_id]               <- Poisson_commit_value_fn2(param = param_9b, x = 0.025, x_star = out_commit[i,'x_star',m_id], x_bounds = temp)
    out_commit[i,paste0('x_',1:n_x),m_id] <- Poisson_commit_value_fn2(param = param_9b, x = x_eval2[[m_id]], x_star = out_commit[i,'x_star',m_id], x_bounds = temp)
    out_commit[i,'dt',m_id]               <- 1/lambda_vec2[i]
    
    # store results
    write.table(out_commit[i,c('dt','m','v0','x_bar','x_bar_bar','x_star'),m_id], file = paste0("poisson_model_commit_maturity_",m_id,"_lambda_",i,".csv"), eol = "\n", na = "NA", fileEncoding = "", col.names = F)
  }
  
  # fill commitment case --- case lambda = +infty
  out_commit[n_lambda2,'x_bar',m_id]     <- (param_9b['r']+param_9b['m'])*(1-param_9b['alpha'])/(param_9b['kappa']+param_9b['m'])/(param_9b['r']-param_9b['mu'])
  out_commit[n_lambda2,'x_bar_bar',m_id] <- out_commit[n_lambda2,'x_bar',m_id]
  out_commit[n_lambda2,'x_star',m_id]    <- out_commit[n_lambda2,'x_bar',m_id]
  out_commit[n_lambda2,'v0',m_id]        <- (out_commit[n_lambda2,'x_bar',m_id])*(param_9b['kappa']+param_9b['m'])/(param_9b['r']+param_9b['m']) + param_9b['alpha']/(param_9b['delta']-param_9b['mu'])
  out_commit[n_lambda2,paste0('x_',1:n_x),m_id] <- (out_commit[n_lambda2,'x_bar',m_id] - x_eval2[[m_id]])*(param_9b['kappa']+param_9b['m'])/(param_9b['r']+param_9b['m']) + param_9b['alpha']/(param_9b['delta']-param_9b['mu'])
  out_commit[n_lambda2,'dt',m_id]        <- 0
  out_commit[n_lambda2,'m',m_id]         <- m_vec[m_id]
}

# export main results
setwd(analysis_dir)
for (m_id in (1:n_m)){
  write.table(out_commit[,c('dt','v0','x_bar','x_bar_bar','x_star'),m_id], 
              file = paste0('commitment_results_m_id_',m_id,'.csv'), eol = "\n", na = "NA", fileEncoding = "", sep = "|", row.names = FALSE)
}

#======================================
# Poisson Issuance Model -- Paper plots
#======================================

# set print directory
setwd(analysis_dir)

# select a particular maturity
m_id          <- 2
param_9['m']  <- m_vec[m_id]

# plot welfare vs. dt for no commitment model
y_min <- 1/(param_9['delta']-param_9['mu'])
y_max <- max(out_no_commit[,'v0',m_id])
pdf("welfare_vs_dt_transformed_no_commit.pdf")
plot(x = 1/(1+lambda_vec), y = out_no_commit[,'v0',m_id], type='l', xlab = '', 
     ylim = c(y_min, y_max), xlim = c(0, 1), axes = F,
     ylab = "",col="blue", lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n')
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par("usr")
arrows(x0 = 0, y0 = y_min, x1 = .15*(u[2]-u[1]), y1 = u[3] + .3*(u[4]-u[3]), 
       angle = 10, code = 1, col = 'black', lty = 1, lwd = thickness, length = .2)
text(x = .14*(u[2]-u[1]), y = u[3] + .35*(u[4]-u[3]), "continuous", cex = cex_plot, pos = 4)
text(x = .17*(u[2]-u[1]), y = u[3] + .3*(u[4]-u[3]), "trading", cex = cex_plot, pos = 4)
arrows(x0 = 1, y0 = y_min, x1 = 1 - .15*(u[2]-u[1]), y1 = u[3] + .5*(u[4]-u[3]), 
       angle = 10, code = 1, col = 'black', lty = 1, lwd = thickness, length = .2)
text(x = 1 - .25*(u[2]-u[1]), y = u[3] + .55*(u[4]-u[3]), "autarky", cex = cex_plot, pos = 4)
axis(2, col="black", lwd=thickness, cex.axis = cex_plot)
axis(1, col="black", lwd=thickness, cex.axis = cex_plot, at = 1/(1+1/c(0,0.5,1,2,3,5,10,1000)), labels = c(0,0.5,1,2,3,5,10,expression(infinity)))
mtext(1,text=expression('expected trading interval'~Delta),line=2.5,cex = cex_plot)
mtext(2,text=expression(v[Delta](0)),line=2.5,cex = cex_plot)
abline(h = y_min, col = 'red', lwd = thickness, lty = 2)
legend("topleft", legend = c('Poisson model', 'no-commitment'), 
       col = c('blue', 'red'), bty="n", lty = c(1,2), cex = cex_plot, lwd = thickness)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()


# plot welfare vs. dt for commitment model
y_min <- 1/(param_9['delta']-param_9['mu'])
y_max <- 1/(param_9['r']-param_9['mu'])
pdf("welfare_vs_dt_transformed_commit.pdf")
plot(x = 1/(1+lambda_vec2), y = out_commit[,'v0',m_id], type='l', xlab = '', 
     ylim = c(y_min, y_max), xlim = c(0, 1), axes = F,
     ylab = "",col="blue", lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n')
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par("usr")
axis(2, col="black", lwd=thickness, cex.axis = cex_plot)
axis(1, col="black", lwd=thickness, cex.axis = cex_plot, at = 1/(1+1/c(0,0.5,1,2,3,5,10,1000)), labels = c(0,0.5,1,2,3,5,10,expression(infinity)))
mtext(1,text=expression('expected trading interval'~Delta),line=2.5,cex = cex_plot)
mtext(2,text=expression(v[Delta~",c"](0)),line=2.5,cex = cex_plot)
abline(h = y_min, col = 'red', lwd = thickness, lty = 2)
abline(h = y_max, col = 'green', lwd = thickness, lty = 3)
legend("right", legend = c('Poisson commit. model', 'no-commitment', 'full-commitment'), 
       col = c('blue', 'red', 'green'), bty="n", lty = 1:3, cex = cex_plot, lwd = thickness)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()


# plot welfare vs. dt for both no commitment and commitment model
y_min <- 1/(param_9['delta']-param_9['mu'])
y_max <- 1/(param_9['r']-param_9['mu'])
pdf("welfare_vs_dt_transformed.pdf")
plot(x = 1/(1+lambda_vec), y = out_no_commit[,'v0',m_id], type='l', xlab = "", 
     ylim = c(y_min, y_max), xlim = c(0, 1), axes = F,
     ylab = "",col="blue", lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n')
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par("usr")
arrows(x0 = 0, y0 = y_min, x1 = .15*(u[2]-u[1]), y1 = u[3] + .3*(u[4]-u[3]), 
       angle = 10, code = 1, col = 'black', lty = 1, lwd = thickness, length = .2)
text(x = .14*(u[2]-u[1]), y = u[3] + .35*(u[4]-u[3]), "continuous", cex = cex_plot, pos = 4)
text(x = .17*(u[2]-u[1]), y = u[3] + .3*(u[4]-u[3]), "trading", cex = cex_plot, pos = 4)
arrows(x0 = 1, y0 = y_min, x1 = 1 - .15*(u[2]-u[1]), y1 = u[3] + .5*(u[4]-u[3]), 
       angle = 10, code = 1, col = 'black', lty = 1, lwd = thickness, length = .2)
text(x = 1 - .25*(u[2]-u[1]), y = u[3] + .55*(u[4]-u[3]), "autarky", cex = cex_plot, pos = 4)
lines(x = 1/(1+lambda_vec2), y = out_commit[,'v0',m_id], col="purple", lwd = thickness)
axis(2, col="black", lwd=thickness, cex.axis = cex_plot)
axis(1, col="black", lwd=thickness, cex.axis = cex_plot, at = 1/(1+1/c(0,0.5,1,2,3,5,10,1000)), labels = c(0,0.5,1,2,3,5,10,expression(infinity)))
mtext(1,text=expression('expected trading interval'~Delta),line=2.5,cex = cex_plot)
mtext(2,text='sovereign welfare',line=2.5,cex = cex_plot)
abline(h = y_min, col = 'red', lwd = thickness, lty = 2)
abline(h = y_max, col = 'green', lwd = thickness, lty = 3)
legend("topright", legend = c('infrequent trading, no commitment', 'infrequent trading, commitment','autarky','first-best'), 
       col = c('blue', 'purple','red','green'), bty="n", lty = c(1,1,2,3), cex = cex_plot, lwd = thickness)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()


# plot default boundary vs. lambda for no-commitment model
y_min <- no_commit_cutoff_fn(param_9)
y_max <- max(out_no_commit[,'x_bar',m_id])
pdf("x_bar_vs_dt_transformed_no_commit.pdf")
plot(x = 1/(1+lambda_vec), y = out_no_commit[,'x_bar',m_id], type='l', xlab = '', ylim = c(y_min, y_max), axes = F,
     ylab = "",col="blue", lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n')
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
abline(h = y_min, col = 'red', lwd = thickness, lty = 2)
legend("topleft", legend = c('Poisson model', 'no-commitment'), 
       col = c('blue', 'red'), bty="n", lty = c(1,2), cex = cex_plot, lwd = thickness)
axis(2, col="black", lwd=thickness, cex.axis = cex_plot)
axis(1, col="black", lwd=thickness, cex.axis = cex_plot, at = 1/(1+1/c(0,0.5,1,2,3,5,10,1000)), labels = c(0,0.5,1,2,3,5,10,expression(infinity)))
mtext(1,text=expression('expected trading interval'~Delta),line=2.5,cex = cex_plot)
mtext(2,text=expression(bar(x)[Delta]),line=2.5,cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()


# plot default boundaries for commitment and no commitment model
y_min <- min(out_commit[,'x_bar_bar',m_id], out_commit[,'x_bar',m_id], out_no_commit[,'x_bar',m_id])
y_max <- 1/(param_9['r']-param_9['mu'])
pdf("x_bar_vs_dt_transformed.pdf")
plot(x = 1/(1+lambda_vec2), y = out_commit[,'x_bar',m_id], type='l', xlab = '', ylim = c(y_min, y_max), lty = 2, 
     ylab = "",col="pink", axes = F, lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n')
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = 1/(1+lambda_vec2), y = out_commit[,'x_bar_bar',m_id], col="purple", lwd = thickness, lty = 3)
lines(x = 1/(1+lambda_vec), y = out_no_commit[,'x_bar',m_id], col="blue", lwd = thickness, lty = 1)
abline(h = 1/(param_9['r'] - param_9['mu']), col = 'green', lwd = thickness, lty = 3)
legend("bottomleft", legend = c(expression(bar(x)[Delta]),expression(bar(x)[Delta][',c']),expression(bar(bar(x))[Delta][',c'])), 
       col = c('blue','pink','purple'), bty="n", lty = 1:3, cex = cex_plot, lwd = thickness)
axis(2, col="black", lwd=thickness, cex.axis = cex_plot)
axis(1, col="black", lwd=thickness, cex.axis = cex_plot, at = 1/(1+1/c(0,0.5,1,2,3,5,10,1000)), labels = c(0,0.5,1,2,3,5,10,expression(infinity)))
mtext(1,text=expression('expected trading interval'~Delta),line=2.5,cex = cex_plot)
mtext(2,text='default boundaries',line=2.5,cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()


# plot target leverage for commitment and no commitment model
y_min <- min(out_commit[,'x_star',m_id], out_no_commit[,'x_a',m_id], no_commit_attraction_pt_fn(param = param_9))
y_max <- max(out_commit[,'x_star',m_id], out_no_commit[,'x_a',m_id], no_commit_attraction_pt_fn(param = param_9))
pdf("target_x_vs_dt_transformed.pdf")
plot(x = 1/(1+lambda_vec2), y = out_commit[,'x_star',m_id], type='l', xlab = '', ylim = c(y_min, y_max), lty = 2,
     ylab = "",col="purple", axes = F, lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n')
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
u <- par('usr')
lines(x = 1/(1+lambda_vec), y = pmin(no_commit_attraction_pt_fn(param = param_9),out_no_commit[,'x_a',m_id]), col="blue", lwd = thickness, lty = 1)
abline(h = no_commit_attraction_pt_fn(param = param_9), col = 'lightblue', lwd = thickness, lty = 3)
abline(h = 1/(param_9['r'] - param_9['mu']), col = 'green', lwd = thickness, lty = 4)
text(x = u[2] - .04*(u[2]-u[1]), y = no_commit_attraction_pt_fn(param = param_9) + .04*(u[4]-u[3]), expression(x[a]), cex = cex_plot, col = 'black') 
legend("bottomleft", legend = c(expression(x[Delta][',a']),expression(x[Delta][',c']^"*")), 
       col = c('blue','purple'), bty="n", lty = 1:2, cex = cex_plot, lwd = thickness)
axis(2, col="black", lwd=thickness, cex.axis = cex_plot)
axis(1, col="black", lwd=thickness, cex.axis = cex_plot, at = 1/(1+1/c(0,0.5,1,2,3,5,10,1000)), labels = c(0,0.5,1,2,3,5,10,expression(infinity)))
mtext(1,text=expression('expected trading interval'~Delta),line=2.5,cex = cex_plot)
mtext(2,text='target debt-to-income',line=2.5,cex = cex_plot)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()


# plot welfare at zero for no commitment and commitment model for different values of m
col_vec<- colours()[(1:n_m)*8+20]
y_min  <- 1/(param_9['delta']-param_9['mu'])
y_max  <- max(max(out_no_commit[,'v0',]),max(out_commit[,'v0',]))
pdf("welfare_vs_dt_transformed_varying_m.pdf")
plot(x = 1/(1+lambda_vec), y = out_no_commit[,'v0',1], xlab = '', 
     ylim = c(y_min, y_max), xlim = c(0, 1), axes = F,
     ylab = "",col=col_vec[1], lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n', 
     type = 'b', pch = 17, lty = 1)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = 1/(1+lambda_vec2), y = out_commit[,'v0',1], lwd = thickness, col = col_vec[1], lty = 1)
u <- par("usr")
for (m_id in (2:n_m)){
  lines(x = 1/(1+lambda_vec), y = out_no_commit[,'v0',m_id], col = col_vec[m_id], lwd = thickness, lty = m_id, pch = 17, type = 'b')
  lines(x = 1/(1+lambda_vec2), y = out_commit[,'v0',m_id], col = col_vec[m_id], lwd = thickness, lty = m_id)
}
axis(2, col="black", lwd=thickness, cex.axis = cex_plot)
axis(1, col="black", lwd=thickness, cex.axis = cex_plot, at = 1/(1+1/c(0,0.5,1,2,3,5,10,1000)), labels = c(0,0.5,1,2,3,5,10,expression(infinity)))
mtext(1,text=expression('expected trading interval'~Delta),line=2.5,cex = cex_plot)
mtext(2,text='sovereign welfare',line=2.5,cex = cex_plot)
abline(h = y_min, col = 'red', lwd = thickness, lty = 5)
abline(h = y_max, col = 'green', lwd = thickness, lty = 6)
legend("topright", legend = paste0('1/m = ',round(1/m_vec,1),' years'), 
       col = col_vec, bty="n", lty = 1:n_m, cex = cex_plot, lwd = thickness)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()


# plot welfare at zero for no commitment and commitment model for different values of m
# this is the .eps version
setEPS()
postscript("welfare_vs_dt_transformed_varying_m.eps")
plot(x = 1/(1+lambda_vec), y = out_no_commit[,'v0',1], xlab = '', 
     ylim = c(y_min, y_max), xlim = c(0, 1), axes = F,
     ylab = "",col=col_vec[1], lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n', 
     type = 'b', pch = 17, lty = 1)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = 1/(1+lambda_vec2), y = out_commit[,'v0',1], lwd = thickness, col = col_vec[1], lty = 1)
u <- par("usr")
for (m_id in (2:n_m)){
  lines(x = 1/(1+lambda_vec), y = out_no_commit[,'v0',m_id], col = col_vec[m_id], lwd = thickness, lty = m_id, pch = 17, type = 'b')
  lines(x = 1/(1+lambda_vec2), y = out_commit[,'v0',m_id], col = col_vec[m_id], lwd = thickness, lty = m_id)
}
axis(2, col="black", lwd=thickness, cex.axis = cex_plot)
axis(1, col="black", lwd=thickness, cex.axis = cex_plot, at = 1/(1+1/c(0,0.5,1,2,3,5,10,1000)), labels = c(0,0.5,1,2,3,5,10,expression(infinity)))
mtext(1,text=expression('expected trading interval'~Delta),line=2.5,cex = cex_plot)
mtext(2,text='sovereign welfare',line=2.5,cex = cex_plot)
abline(h = y_min, col = 'red', lwd = thickness, lty = 5)
abline(h = y_max, col = 'green', lwd = thickness, lty = 6)
legend("topright", legend = paste0('1/m = ',round(1/m_vec,1),' years'), 
       col = col_vec, bty="n", lty = 1:n_m, cex = cex_plot, lwd = thickness)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()


# plot welfare at zero for no commitment and commitment model for presentation
y_min  <- 1/(param_9['delta']-param_9['mu'])
y_max  <- max(max(out_no_commit[,'v0',]),max(out_commit[,'v0',]))
pdf("welfare_vs_dt_transformed_presentation.pdf")
plot(x = 1/(1+lambda_vec), y = out_no_commit[,'v0',2], type='l', xlab = '', 
     ylim = c(y_min, y_max), xlim = c(0, 1), axes = F,
     ylab = "",col='blue', lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, bty = 'n')
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = 1/(1+lambda_vec2), y = out_commit[,'v0',2], lwd = thickness, lty = 2, col = 'purple')
u <- par("usr")
axis(2, col="black", lwd=thickness, cex.axis = cex_plot)
axis(1, col="black", lwd=thickness, cex.axis = cex_plot, at = 1/(1+1/c(0,0.5,1,2,3,5,10,1000)), labels = c(0,0.5,1,2,3,5,10,expression(infinity)))
mtext(1,text=expression('expected trading interval'~Delta),line=2.5,cex = cex_plot)
mtext(2,text='sovereign welfare',line=2.5,cex = cex_plot)
abline(h = y_min, col = 'red', lwd = thickness, lty = 3)
abline(h = y_max, col = 'green', lwd = thickness, lty = 3)
arrows(x0 = u[1] + .2*(u[2]-u[1]), y0 = y_min, x1 = u[1] + .2*(u[2]-u[1]), y1 = u[3] + .3*(u[4]-u[3]), 
       angle = 10, code = 1, col = 'black', lty = 1, lwd = thickness, length = .15)
text(x = u[1] + .1*(u[2]-u[1]), y = u[3] + .4*(u[4]-u[3]), "autarky", cex = cex_plot, pos = 4)
text(x = u[1] + .1*(u[2]-u[1]), y = u[3] + .35*(u[4]-u[3]), "welfare", cex = cex_plot, pos = 4)
arrows(x0 = u[1] + .4*(u[2]-u[1]), y0 = y_max, x1 = u[1] + .4*(u[2]-u[1]), y1 = u[3] + .8*(u[4]-u[3]), 
       angle = 10, code = 1, col = 'black', lty = 1, lwd = thickness, length = .15)
text(x = u[1] + .3*(u[2]-u[1]), y = u[3] + .75*(u[4]-u[3]), "first best", cex = cex_plot, pos = 4)
legend("right", legend = c('no commitment','commitment'), 
       col = c('blue', 'purple'), bty="n", lty = c(1,2), cex = cex_plot, lwd = thickness)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

#=======================================
# Compute equilibrium with risk-aversion
#=======================================

# create output directory for pictures for illustrative plots
setwd(output_dir)
if (!file.exists('risk_aversion')) dir.create(file.path(output_dir, 'risk_aversion'))

# store output in appropriate folder
analysis_dir    <- file.path(output_dir, 'risk_aversion')
setwd(analysis_dir)

# economic parameters
param_risk_av              <- c()
param_risk_av['mu']        <- .02     # income growth rate
param_risk_av['sigma']     <- .1      # income volatility
param_risk_av['m']         <- 1/20    # bond maturity intensity
param_risk_av['theta']     <- .5      # debt haircut upon default
param_risk_av['alpha']     <- .96     # GDP drop upon default
param_risk_av['r']         <- .05     # international risk free rates
param_risk_av['kappa']     <- param_0['r']      # coupon rate
param_risk_av['delta']     <- param_0['r']+.02  # impatience of government
param_risk_av['nu']        <- 0       # international risk prices
param_risk_av['gamma']     <- 2       # risk aversion
param_risk_av['eta']       <- 0       # tax on issuance proceeds

# numerical parameters
num_param_risk_av               <- c()      # declare numerical parameters
num_param_risk_av['n_x']        <- 2000     # number of discretization points
num_param_risk_av['dt']         <- .25      # time interval for false transient method
num_param_risk_av['max_iter']   <- 1000     # max number of iterations
num_param_risk_av['tol']        <- .00001   # tolerance
num_param_risk_av['relax']      <- .02      # relaxation parameter: policy = relax*new policy+(1-relax)*old 
num_param_risk_av['x_max']      <- 2        # max debt-to-income
num_param_risk_av['print_flag'] <- F        # print flag for error checking purposes
num_param_risk_av['policy']     <- 2        # what policy to use? choice: 0 = optimal, 1 = risk-neutral, 2 = no issuances

# plotting parameters
plot_risk_av                     <- c()
plot_risk_av[['n']]              <- 200
plot_risk_av[['scen_name']]      <- 'base_case'
plot_risk_av[['output_dir']]     <- analysis_dir
plot_risk_av[['line_thickness']] <- thickness
plot_risk_av[['font_size']]      <- cex_plot
plot_risk_av[['use_num_sol']]    <- T

#========================================
# Compute comparative static w.r.t. gamma
#========================================

# create vector of risk aversion, as well as numerical parameters for each gamma parameter
gamma_vec    <- c(.025,.05,.1,.25,.5,.75,.95,1.05,1.25,1.5,1.75,2)
tol_vec      <- c(.000001,.00001,.00001,.00001,.00001,.00001,.00001,.00001,.00001,.00001,.00001,.00001)
n_gamma      <- length(gamma_vec)

# pre-compute risk-neutral policies
RN_policy <- list()
RN_policy[['x_bar']] <- no_commit_cutoff_fn(param = param_risk_av)
RN_policy[['x']]     <- seq(from = 0.2, to = RN_policy[['x_bar']], length.out = 200)
RN_policy[['g']]     <- no_commit_g_fn(param = param_risk_av, x = RN_policy[['x']])

# parameters when gamma = 0
param_RN          <- param_risk_av
param_RN['gamma'] <- 0

# launch computations
for (i in (1:n_gamma)){
  print(paste0('percentage completed = ',100*(i-1)/n_gamma,'%'))
  param_risk_av['gamma']          <- gamma_vec[i]
  num_param_risk_av['tol']        <- tol_vec[i]
  
  # create folder where to store optimal behavior for gamma > 0
  setwd(analysis_dir)
  if (!file.exists(as.character(gamma_vec[i]))) dir.create(file.path(analysis_dir, as.character(gamma_vec[i])))
  setwd(file.path(analysis_dir, as.character(gamma_vec[i])))
  plot_risk_av[['output_dir']] <- file.path(analysis_dir, as.character(gamma_vec[i]))
  
  # solve for equilibrium for gamma > 0
  sol_risk_av <- risk_aversion_fn(param = param_risk_av, num_param = num_param_risk_av, policies = NA)
  eq_plot_fn(param = param_risk_av, eq_sol = sol_risk_av, plot_param = plot_risk_av, num_param = num_param_risk_av)
  
  # store optimal policy for risk-averse borrower
  risk_av_policy            <- list()
  risk_av_policy[['x_bar']] <- sol_risk_av$other_moments['x_bar']
  risk_av_policy[['x']]     <- sol_risk_av$eq_objects$x
  risk_av_policy[['g']]     <- sol_risk_av$eq_objects$g
  
  # create folder where to store suboptimal behavior for gamma > 0 (where policies are RN policies)
  setwd(analysis_dir)
  if (!file.exists(as.character(paste0(gamma_vec[i],' - RN policies')))) dir.create(file.path(analysis_dir, as.character(paste0(gamma_vec[i],' - RN policies'))))
  setwd(file.path(analysis_dir, as.character(paste0(gamma_vec[i],' - RN policies'))))
  plot_risk_av[['output_dir']] <- file.path(analysis_dir, as.character(paste0(gamma_vec[i],' - RN policies')))
  
  # solve for equilibrium for gamma > 0 where policies are RN policies
  sol_risk_av <- risk_aversion_fn(param = param_risk_av, num_param = num_param_risk_av, policies = RN_policy)
  eq_plot_fn(param = param_risk_av, eq_sol = sol_risk_av, plot_param = plot_risk_av, num_param = num_param_risk_av)
  
  # create folder where to store suboptimal behavior for gamma = 0
  setwd(analysis_dir)
  if (!file.exists(as.character(paste0(gamma_vec[i],' - RN borrower')))) dir.create(file.path(analysis_dir, as.character(paste0(gamma_vec[i],' - RN borrower'))))
  setwd(file.path(analysis_dir, as.character(paste0(gamma_vec[i],' - RN borrower'))))
  plot_risk_av[['output_dir']] <- file.path(analysis_dir, as.character(paste0(gamma_vec[i],' - RN borrower')))
  
  # solve for equilibrium for gamma > 0 where policy is zero debt issuance
  sol_risk_av <- risk_aversion_fn(param = param_RN, num_param = num_param_risk_av, policies = risk_av_policy)
  eq_plot_fn(param = param_RN, eq_sol = sol_risk_av, plot_param = plot_risk_av, num_param = num_param_risk_av)
}

#=============================================================
# Plots for risk aversion sensitivity -- plot entire functions
#=============================================================

# collect results
setwd(analysis_dir)

# store key outputs
gamma_plot_vec   <- c(.05,.1,.25,.5)
n_gamma_plot     <- length(gamma_plot_vec)
out_names        <- c("x_bar","default_rate","d_recov","cons_mult","v_aut","v_0")
x_bar_vec        <- array(NA,n_gamma_plot)

# store functions
out_fn           <- list()
for (i in (1:n_gamma_plot)){
  setwd(file.path(analysis_dir, as.character(gamma_plot_vec[i])))
  # download functions
  out_fn[[i]]      <- list()
  prout            <- read.csv("policies.csv", sep = " ", header = T)
  for (field in names(prout)){
    out_fn[[i]][[field]] <- prout[,field]
  }
  # download default boundary
  prout            <- read.csv("other_moments.csv", sep = " ", header = F)
  rownames(prout)  <- prout[,1]
  prout            <- prout %>% dplyr::select(-1)
  x_bar_vec[i]     <- prout['x_bar',1]
}

# store risk-neutral case separately -- functions
setwd(analysis_dir)
if (!file.exists(as.character(0))) dir.create(file.path(analysis_dir, as.character(0)))
setwd(file.path(analysis_dir, as.character(0)))
RN_sol            <- list()
x_bar_RN          <- no_commit_cutoff_fn(param_risk_av)
RN_sol[['x']]     <- seq(from = .01, to = x_bar_RN, length = 1000)
RN_sol[['v']]     <- no_commit_value_fn(param_risk_av, x = RN_sol[['x']])
RN_sol[['d']]     <- no_commit_debt_px_fn(param_risk_av, x = RN_sol[['x']])
RN_sol[['g']]     <- no_commit_g_fn(param_risk_av, x = RN_sol[['x']])
RN_sol[['cons']]  <- no_commit_c_fn(param_risk_av, x = RN_sol[['x']])
RN_sol[['spread']]<- no_commit_credit_spread_fn(param_risk_av, x = RN_sol[['x']])
RN_sol[['dist']]  <- no_commit_ergodic_dist_fn(param_risk_av)[["density"]](RN_sol[['x']])

# plots of policy and value functions
setwd(analysis_dir)
fn_to_plot    <- c('cons','d','spread','g','dist','v')
title_plot    <- c('consumption-to-income c(x)', 'debt price d(x)', 'credit spread', 'issuance rate g(x)', 'ergodic density', '(scaled) value function v(x)')
names(title_plot) <- fn_to_plot
col_vec       <- colours()[20+(1:n_gamma_plot)*7]
y_min         <- c(.8,.4,0,0,0,15)
y_max         <- c(1.5,1,.2,.8,8,35)
names(y_min)  <- fn_to_plot
names(y_max)  <- fn_to_plot
leg_pos       <- c('bottomleft', 'bottomleft', 'topleft', 'bottomleft', 'topleft', 'topleft')
names(leg_pos)<- fn_to_plot
for (fn_plot in fn_to_plot){
  # print both in pdf and eps version
  for (file_format in c('pdf','eps')){
    if (file_format == 'pdf') {
      lty_count <- 2
      pdf(paste0('risk_aversion_',fn_plot,'.pdf'))
    } else {
      lty_count <- 2
      setEPS()
      postscript(paste0('risk_aversion_',fn_plot,'.eps'))
    }
    par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
    plot(x = RN_sol[['x']], y = RN_sol[[fn_plot]], bty = 'n', type = 'l',axes = 'F', main = '', ylim = c(y_min[fn_plot],y_max[fn_plot]), 
         xlab = '', ylab = '', col = 'black', lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, lty = 1)
    for (gam in (1:n_gamma_plot)){
      lines(x = out_fn[[gam]][['x']][out_fn[[gam]][['x']]<.995*x_bar_vec[gam]], 
            y = out_fn[[gam]][[fn_plot]][out_fn[[gam]][['x']]<.995*x_bar_vec[gam]], 
            col = col_vec[gam], lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, lty = lty_count)
      lty_count <- lty_count + 1
    }
    grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
    axis(2, col="black", lwd = 1, cex.axis = cex_plot)
    mtext(2,text = title_plot[fn_plot], line = 2, cex = cex_plot)
    axis(1, col="black", lwd = 1, cex.axis = cex_plot)
    mtext(1,text="debt-to-income",line=2, cex = cex_plot)
    legend(leg_pos[fn_plot], col = c('black',col_vec), bty="n", lty = (1:(n_gamma_plot+1)), cex = cex_plot, lwd = thickness, 
           legend = sapply(1:(n_gamma_plot+1), function(x) as.expression(substitute(gamma == B,
                                                                                    list(B = c(0,gamma_plot_vec)[x])))))
    box(which = "plot", bty = "l", lwd = thickness)
    dev.off() 
  }
}

#=======================================================================
# Plots for risk aversion sensitivity -- plot moments and other outcomes
#=======================================================================

# collect results
setwd(analysis_dir)

# store key outputs
gamma_plot_vec   <- c(.025,.05,.1,.25,.5,.75,.95,1.05,1.25,1.5,1.75,2)
n_gamma_plot     <- length(gamma_plot_vec)
output           <- array(NA, dim = c(n_gamma_plot,11))
out_names        <- c("x_bar","default_rate","d_recov","cons_mult","v_aut","v_0")
colnames(output) <- c("gamma",out_names,"v_0_RN_policies","cons_mult_RN_policies","v_0_RN_borrower","cons_mult_RN_borrower")
out_legend       <- c("default boundary",
                      "default rate",
                      "debt recovery rate",
                      "consumption multiple",
                      "(scaled) autarky value",
                      "(scaled) welfare")
names(out_legend)<- out_names
moments          <- array(NA, dim = c(n_gamma_plot,7))
moments_names    <- c("x_mean","x_stdev","v_mean","v_stdev","spread_mean","spread_stdev")
colnames(moments)<- c("gamma",moments_names)
mom_legend       <- c("debt-to-income average",
                      "debt-to-income standard deviation",
                      "value function average",
                      "value function standard deviation",
                      "credit spread average",
                      "credit spread standard deviation")
names(mom_legend) <- moments_names

# store outcomes
out_fn <- list()
for (i in (1:n_gamma_plot)){
  setwd(file.path(analysis_dir, gamma_plot_vec[i]))
  # download other moments
  prout            <- read.csv("other_moments.csv", sep = " ", header = F)
  rownames(prout)  <- prout[,1]
  prout            <- prout %>% dplyr::select(-1)
  output[i,"gamma"]<- gamma_plot_vec[i]
  for (out_name in out_names){
    output[i,out_name] <- prout[out_name,1]
  }
  
  # download moments
  prout                    <- read.csv("moments.csv", sep = " ", header = T)
  moments[i,"gamma"]       <- gamma_plot_vec[i]
  moments[i,"x_mean"]      <- prout['x','mean']
  moments[i,"x_stdev"]     <- prout['x','stdev']
  moments[i,"v_mean"]      <- prout['v','mean']
  moments[i,"v_stdev"]     <- prout['v','stdev']
  moments[i,"spread_mean"] <- prout['spread','mean']
  moments[i,"spread_stdev"]<- prout['spread','stdev']
  
  # download v0 for RN policies
  setwd(file.path(analysis_dir, paste0(gamma_plot_vec[i]," - RN policies")))
  prout            <- read.csv("other_moments.csv", sep = " ", header = F)
  rownames(prout)  <- prout[,1]
  prout            <- prout %>% dplyr::select(-1)
  output[i,"v_0_RN_policies"]      <- prout["v_0",1]
  output[i,"cons_mult_RN_policies"]<- (output[i,"v_0"]/output[i,"v_0_RN_policies"])^(1/(1-gamma_plot_vec[i]))
  
  # download v0 for RN borrower and risk-averse policies
  setwd(file.path(analysis_dir, paste0(gamma_plot_vec[i]," - RN borrower")))
  prout            <- read.csv("other_moments.csv", sep = " ", header = F)
  rownames(prout)  <- prout[,1]
  prout            <- prout %>% dplyr::select(-1)
  output[i,"v_0_RN_borrower"]      <- prout["v_0",1]
  output[i,"cons_mult_RN_borrower"]<- output[i,"v_0_RN_borrower"]*(param_risk_av['delta'] - param_risk_av['mu'])
}

# store risk-neutral case separately -- store values
setwd(analysis_dir)
if (!file.exists(as.character(0))) dir.create(file.path(analysis_dir, as.character(0)))
setwd(file.path(analysis_dir, as.character(0)))
RN_out            <- array(NA,6)
names(RN_out)     <- out_names
RN_out['x_bar']   <- no_commit_cutoff_fn(param_risk_av)
RN_out['default_rate'] <- no_commit_ergodic_dist_fn(param_risk_av)[["default_rate"]]
RN_out['d_recov'] <- no_commit_debt_px_fn(param_risk_av, x = RN_out['x_bar'])
RN_out['cons_mult'] <- 1
RN_out['v_aut']   <- 1/(param_risk_av['delta'] - param_risk_av['mu'])
RN_out['v_0']     <- RN_out['v_aut']

# plots of moments
setwd(analysis_dir)
for (mom in moments_names){
  pdf(paste0('risk_aversion_',mom,'.pdf'))
  par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
  plot(x = gamma_plot_vec, y = moments[,mom], bty = 'n', type = 'l',axes = 'F', main = '', 
       xlab = '', ylab = '', col = 'black', lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  axis(2, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(2,text = mom_legend[mom], line = 2, cex = cex_plot)
  axis(1, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(1,text=expression(gamma),line=2, cex = cex_plot)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
}

# plots of other outcomes
setwd(analysis_dir)
for (outcome in out_names){
  pdf(paste0('risk_aversion_',outcome,'.pdf'))
  par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
  plot(x = c(0,gamma_plot_vec), y = c(RN_out[outcome], output[,outcome]), bty = 'n', type = 'l',axes = 'F', main = '', 
       xlab = '', ylab = '', col = 'black', lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  axis(2, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(2,text = out_legend[outcome], line = 2, cex = cex_plot)
  axis(1, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(1,text=expression(gamma),line=2, cex = cex_plot)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
}

# plot v_0 vs. v_0 suboptimal (using RN policies)
pdf('risk_aversion_v_0_suboptimal.pdf')
par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
plot(x = c(0,gamma_plot_vec[gamma_plot_vec<1]), y = c(RN_out['v_0'], output[gamma_plot_vec<1,'v_0']), 
     bty = 'n', type = 'l',axes = 'F', main = '', 
     xlab = '', ylab = '', col = 'black', lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, 
     ylim = c(min(output[,'v_0']),max(output[,'v_0'])), xlim = c(0,max(gamma_plot_vec)))
lines(x = gamma_plot_vec[gamma_plot_vec>=1], y = output[gamma_plot_vec>=1,'v_0'], col = 'black', lwd = thickness)
lines(x = c(0,gamma_plot_vec[gamma_plot_vec<1]), y = c(RN_out['v_0'], output[gamma_plot_vec<1,'v_0_RN_policies']), 
      col = 'red', lwd = thickness, lty = 2)
lines(x = gamma_plot_vec[gamma_plot_vec>=1], y = output[gamma_plot_vec>=1,'v_0_RN_policies'], 
      col = 'red', lwd = thickness, lty = 2)
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
axis(2, col="black", lwd = 1, cex.axis = cex_plot)
mtext(2,text = "v(0)", line = 2, cex = cex_plot)
axis(1, col="black", lwd = 1, cex.axis = cex_plot)
mtext(1,text=expression(gamma),line=2, cex = cex_plot)
par(new=T)
plot(x = c(0,gamma_plot_vec), y = c(1,output[,'cons_mult_RN_policies']), axes=F, xlab="", ylab="", 
     type="l", main="", col='red', lty=3, lwd = thickness)
axis(4, col="black", lwd = 1, cex.axis = cex_plot)
mtext(4,text = "consumption multiple", line = 2, cex = cex_plot)
legend("topleft", col = c('black','red','red'), bty="n", lty = c(1,2,3), cex = cex_plot,
       legend = c('v(0)', 'hat_v(0)','consumption multiple'), lwd = thickness)
dev.off()

# plot v_0 with RN borrowers as a function of risk-averse policies
pdf('risk_aversion_v_0_RN_borrower_risk_av_policies.pdf')
par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
plot(x = c(0,gamma_plot_vec), y = c(1, output[,'cons_mult_RN_borrower']), bty = 'n', type = 'l',axes = 'F', main = '', 
     xlab = '', ylab = '', col = 'purple', lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, lty = 2,
     ylim = c(min(c(1,output[,'cons_mult_RN_borrower'],output[,'cons_mult'])),max(c(1,output[,'cons_mult_RN_borrower'],output[,'cons_mult']))), xlim = c(0,max(gamma_plot_vec)))
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = c(0,gamma_plot_vec), y = c(RN_out['cons_mult'], output[,'cons_mult']), col = 'blue', lwd = thickness, lty = 1)
axis(2, col="black", lwd = 1, cex.axis = cex_plot)
mtext(2,text = "consumption multiple", line = 2, cex = cex_plot)
axis(1, col="black", lwd = 1, cex.axis = cex_plot)
mtext(1,text=expression(gamma),line=2, cex = cex_plot)
legend("bottomright", col = c('blue','purple'), bty="n", lty = 1:2, cex = cex_plot,
       legend = c(expression(k[gamma]),expression(hat(k)[gamma])), lwd = thickness)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

# plot v_0 with RN borrowers as a function of risk-averse policies
# plot eps version for paper
setEPS()
postscript('risk_aversion_v_0_RN_borrower_risk_av_policies.eps')
par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
plot(x = c(0,gamma_plot_vec), y = c(1, output[,'cons_mult_RN_borrower']), bty = 'n', type = 'l',axes = 'F', main = '', 
     xlab = '', ylab = '', col = 'purple', lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot, lty = 2,
     ylim = c(min(c(1,output[,'cons_mult_RN_borrower'],output[,'cons_mult'])),max(c(1,output[,'cons_mult_RN_borrower'],output[,'cons_mult']))), xlim = c(0,max(gamma_plot_vec)))
grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
lines(x = c(0,gamma_plot_vec), y = c(RN_out['cons_mult'], output[,'cons_mult']), col = 'blue', lwd = thickness, lty = 1)
axis(2, col="black", lwd = 1, cex.axis = cex_plot)
mtext(2,text = "consumption multiple", line = 2, cex = cex_plot)
axis(1, col="black", lwd = 1, cex.axis = cex_plot)
mtext(1,text=expression(gamma),line=2, cex = cex_plot)
legend("bottomright", col = c('blue','purple'), bty="n", lty = 1:2, cex = cex_plot,
       legend = c(expression(k[gamma]),expression(hat(k)[gamma])), lwd = thickness)
box(which = "plot", bty = "l", lwd = thickness)
dev.off()

#====================================
# Compute comparative static w.r.t. m
#====================================

# create output directory for pictures for illustrative plots
setwd(output_dir)
if (!file.exists('risk_aversion_vs_m')) dir.create(file.path(output_dir, 'risk_aversion_vs_m'))

# store output in appropriate folder
analysis_dir    <- file.path(output_dir, 'risk_aversion_vs_m')
setwd(analysis_dir)

# create vector of risk aversion, as well as numerical parameters for each gamma parameter
m_vec    <- seq(from=.025,to=.1,by=.005)
n_m      <- length(m_vec)
tol_vec  <- rep(.00001,n_m)

# launch computations
for (i in (1:n_m)){
  print(paste0('percentage completed = ',100*(i-1)/n_m,'%'))
  param_risk_av['m']          <- m_vec[i]
  num_param_risk_av['tol']    <- tol_vec[i]
  
  # create folder where to store optimal behavior
  setwd(analysis_dir)
  if (!file.exists(as.character(m_vec[i]))) dir.create(file.path(analysis_dir, as.character(m_vec[i])))
  setwd(file.path(analysis_dir, as.character(m_vec[i])))
  plot_risk_av[['output_dir']] <- file.path(analysis_dir, as.character(m_vec[i]))
  
  # solve for equilibrium for gamma > 0
  sol_risk_av <- risk_aversion_fn(param = param_risk_av, num_param = num_param_risk_av, policies = NA)
  eq_plot_fn(param = param_risk_av, eq_sol = sol_risk_av, plot_param = plot_risk_av, num_param = num_param_risk_av)
}

#==================================================================
# Plots for maturity sensitivity -- plot moments and other outcomes
#==================================================================

# collect results
setwd(analysis_dir)

# store key outputs
m_plot_vec       <- m_vec
n_m_plot         <- length(m_plot_vec)
output           <- array(NA, dim = c(n_m_plot,7))
out_names        <- c("x_bar","default_rate","d_recov","cons_mult","v_aut","v_0")
colnames(output) <- c("m",out_names)
out_legend       <- c("default boundary",
                      "default rate",
                      "debt recovery rate",
                      "consumption multiple",
                      "(scaled) autarky value",
                      "(scaled) welfare")
names(out_legend)<- out_names
moments          <- array(NA, dim = c(n_m,7))
moments_names    <- c("x_mean","x_stdev","v_mean","v_stdev","spread_mean","spread_stdev")
colnames(moments)<- c("m",moments_names)
mom_legend       <- c("debt-to-income average",
                      "debt-to-income standard deviation",
                      "value function average",
                      "value function standard deviation",
                      "credit spread average",
                      "credit spread standard deviation")
names(mom_legend) <- moments_names

# store outcomes
out_fn <- list()
for (i in (1:n_m_plot)){
  setwd(file.path(analysis_dir, m_plot_vec[i]))
  # download other moments
  prout            <- read.csv("other_moments.csv", sep = " ", header = F)
  rownames(prout)  <- prout[,1]
  prout            <- prout %>% dplyr::select(-1)
  output[i,"m"]    <- m_plot_vec[i]
  for (out_name in out_names){
    output[i,out_name] <- prout[out_name,1]
  }
  
  # download moments
  prout                    <- read.csv("moments.csv", sep = " ", header = T)
  moments[i,"m"]           <- m_plot_vec[i]
  moments[i,"x_mean"]      <- prout['x','mean']
  moments[i,"x_stdev"]     <- prout['x','stdev']
  moments[i,"v_mean"]      <- prout['v','mean']
  moments[i,"v_stdev"]     <- prout['v','stdev']
  moments[i,"spread_mean"] <- prout['spread','mean']
  moments[i,"spread_stdev"]<- prout['spread','stdev']
}

# plots of moments
setwd(analysis_dir)
for (mom in moments_names){
  pdf(paste0('risk_aversion_sensi_m_',mom,'.pdf'))
  par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
  plot(x = 1/m_plot_vec, y = moments[,mom], bty = 'n', type = 'l',axes = 'F', main = '', 
       xlab = '', ylab = '', col = 'black', lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  axis(2, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(2,text = mom_legend[mom], line = 2, cex = cex_plot)
  axis(1, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(1,text='debt average life 1/m (years)',line=2, cex = cex_plot)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
}

# plots of other outcomes
setwd(analysis_dir)
for (outcome in out_names){
  pdf(paste0('risk_aversion_sensi_m_',outcome,'.pdf'))
  par(mar=c(4, 4, 4, 4), mgp=c(3, 0.75, 0))
  plot(x = 1/m_plot_vec, y = output[,outcome], bty = 'n', type = 'l',axes = 'F', main = '', 
       xlab = '', ylab = '', col = 'black', lwd = thickness, cex.lab = cex_plot, cex.axis = cex_plot)
  grid(nx = NULL, ny = NULL, col = "lightgray", lty = 1, lwd = .5*thickness, equilogs = TRUE)
  axis(2, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(2,text = out_legend[outcome], line = 2, cex = cex_plot)
  axis(1, col="black", lwd = 1, cex.axis = cex_plot)
  mtext(1,text='debt average life 1/m (years)',line=2, cex = cex_plot)
  box(which = "plot", bty = "l", lwd = thickness)
  dev.off()
}

#======================================
# Move paper plots into separate folder
#======================================

# create export folder
setwd(output_dir)
if (!file.exists('paper_plots')) dir.create(file.path(output_dir, 'paper_plots'))

# copy files, one at a time
file.copy(from = file.path(output_dir, 'base_case/value.eps'), to = file.path(output_dir, 'paper_plots/value.eps'))
file.copy(from = file.path(output_dir, 'base_case/g.eps'), to = file.path(output_dir, 'paper_plots/g.eps'))
file.copy(from = file.path(output_dir, 'base_case/debt.eps'), to = file.path(output_dir, 'paper_plots/debt.eps'))
file.copy(from = file.path(output_dir, 'base_case/spread.eps'), to = file.path(output_dir, 'paper_plots/spread.eps'))

file.copy(from = file.path(output_dir, 'simulation/debt_and_income.eps'), to = file.path(output_dir, 'paper_plots/debt_and_income.eps'))
file.copy(from = file.path(output_dir, 'simulation/debt_to_income.eps'), to = file.path(output_dir, 'paper_plots/debt_to_income.eps'))

file.copy(from = file.path(output_dir, 'maturity_sensitivity/citizen_welfare_loss_vs_m.eps'), to = file.path(output_dir, 'paper_plots/citizen_welfare_loss_vs_m.eps'))
file.copy(from = file.path(output_dir, 'hat_delta_sensitivity/citizen_welfare_loss_vs_hat_delta.eps'), to = file.path(output_dir, 'paper_plots/citizen_welfare_loss_vs_hat_delta.eps'))

file.copy(from = file.path(output_dir, 'markov_switching/debt_price_markov_switching.eps'), to = file.path(output_dir, 'paper_plots/debt_price_markov_switching.eps'))
file.copy(from = file.path(output_dir, 'markov_switching/g_markov_switching.eps'), to = file.path(output_dir, 'paper_plots/g_markov_switching.eps'))

file.copy(from = file.path(output_dir, 'debt_ceiling/debt_price_reflecting_constrained.eps'), to = file.path(output_dir, 'paper_plots/debt_price_reflecting_constrained.eps'))
file.copy(from = file.path(output_dir, 'debt_ceiling/debt_price_reflecting_constrained2.eps'), to = file.path(output_dir, 'paper_plots/debt_price_reflecting_constrained2.eps'))
file.copy(from = file.path(output_dir, 'debt_ceiling/value_function_reflecting_constrained.eps'), to = file.path(output_dir, 'paper_plots/value_function_reflecting_constrained.eps'))
file.copy(from = file.path(output_dir, 'debt_ceiling/debt_price_smooth_constrained.eps'), to = file.path(output_dir, 'paper_plots/debt_price_smooth_constrained.eps'))
file.copy(from = file.path(output_dir, 'debt_ceiling/welfare_analysis/welfare_change_debt_ceiling_1.eps'), to = file.path(output_dir, 'paper_plots/welfare_change_debt_ceiling_1.eps'))

file.copy(from = file.path(output_dir, 'issuance_cap/consumption_issuance_cap.eps'), to = file.path(output_dir, 'paper_plots/consumption_issuance_cap.eps'))
file.copy(from = file.path(output_dir, 'issuance_cap/debt_price_issuance_cap.eps'), to = file.path(output_dir, 'paper_plots/debt_price_issuance_cap.eps'))
file.copy(from = file.path(output_dir, 'issuance_cap/g_issuance_cap.eps'), to = file.path(output_dir, 'paper_plots/g_issuance_cap.eps'))
file.copy(from = file.path(output_dir, 'issuance_cap/value_function_issuance_cap.eps'), to = file.path(output_dir, 'paper_plots/value_function_issuance_cap.eps'))
file.copy(from = file.path(output_dir, 'issuance_cap/welfare_change_issuance_cap_1.eps'), to = file.path(output_dir, 'paper_plots/welfare_change_issuance_cap_1.eps'))

file.copy(from = file.path(output_dir, 'poisson_model/welfare_vs_dt_transformed_varying_m.eps'), to = file.path(output_dir, 'paper_plots/welfare_vs_dt_transformed_varying_m.eps'))

file.copy(from = file.path(output_dir, 'risk_aversion/risk_aversion_v.eps'), to = file.path(output_dir, 'paper_plots/risk_aversion_v.eps'))
file.copy(from = file.path(output_dir, 'risk_aversion/risk_aversion_g.eps'), to = file.path(output_dir, 'paper_plots/risk_aversion_g.eps'))
file.copy(from = file.path(output_dir, 'risk_aversion/risk_aversion_v_0_RN_borrower_risk_av_policies.eps'), to = file.path(output_dir, 'paper_plots/risk_aversion_v_0_RN_borrower_risk_av_policies.eps'))

#=========================================
# Move appendix plots into separate folder
#=========================================

# create export folder
setwd(output_dir)
if (!file.exists('appendix_plots')) dir.create(file.path(output_dir, 'appendix_plots'))

# copy files, one at a time
file.copy(from = file.path(output_dir, 'plots_for_proof/no_debt_buyback.pdf'), to = file.path(output_dir, 'appendix_plots/no_debt_buyback.pdf'))
file.copy(from = file.path(output_dir, 'plots_for_proof/no_debt_issuance.pdf'), to = file.path(output_dir, 'appendix_plots/no_debt_issuance.pdf'))

file.copy(from = file.path(output_dir, 'comparative_statics/debt_to_income_sensi_m.pdf'), to = file.path(output_dir, 'appendix_plots/debt_to_income_sensi_m.pdf'))
file.copy(from = file.path(output_dir, 'comparative_statics/debt_to_income_sensi_r.pdf'), to = file.path(output_dir, 'appendix_plots/debt_to_income_sensi_r.pdf'))
file.copy(from = file.path(output_dir, 'comparative_statics/debt_to_income_sensi_alpha.pdf'), to = file.path(output_dir, 'appendix_plots/debt_to_income_sensi_alpha.pdf'))
file.copy(from = file.path(output_dir, 'comparative_statics/debt_to_income_sensi_theta.pdf'), to = file.path(output_dir, 'appendix_plots/debt_to_income_sensi_theta.pdf'))
file.copy(from = file.path(output_dir, 'comparative_statics/debt_to_income_sensi_mu.pdf'), to = file.path(output_dir, 'appendix_plots/debt_to_income_sensi_mu.pdf'))
file.copy(from = file.path(output_dir, 'comparative_statics/debt_to_income_sensi_sigma.pdf'), to = file.path(output_dir, 'appendix_plots/debt_to_income_sensi_sigma.pdf'))

file.copy(from = file.path(output_dir, 'comparative_statics/deft_rate_and_spreads_sensi_m.pdf'), to = file.path(output_dir, 'appendix_plots/deft_rate_and_spreads_sensi_m.pdf'))
file.copy(from = file.path(output_dir, 'comparative_statics/deft_rate_and_spreads_sensi_r.pdf'), to = file.path(output_dir, 'appendix_plots/deft_rate_and_spreads_sensi_r.pdf'))
file.copy(from = file.path(output_dir, 'comparative_statics/deft_rate_and_spreads_sensi_alpha.pdf'), to = file.path(output_dir, 'appendix_plots/deft_rate_and_spreads_sensi_alpha.pdf'))
file.copy(from = file.path(output_dir, 'comparative_statics/deft_rate_and_spreads_sensi_theta.pdf'), to = file.path(output_dir, 'appendix_plots/deft_rate_and_spreads_sensi_theta.pdf'))
file.copy(from = file.path(output_dir, 'comparative_statics/deft_rate_and_spreads_sensi_mu.pdf'), to = file.path(output_dir, 'appendix_plots/deft_rate_and_spreads_sensi_mu.pdf'))
file.copy(from = file.path(output_dir, 'comparative_statics/deft_rate_and_spreads_sensi_sigma.pdf'), to = file.path(output_dir, 'appendix_plots/deft_rate_and_spreads_sensi_sigma.pdf'))

file.copy(from = file.path(output_dir, 'poisson_model/no_commit_maturity_2/poisson_5/poisson_d.pdf'), to = file.path(output_dir, 'appendix_plots/poisson_d_high.pdf'))
file.copy(from = file.path(output_dir, 'poisson_model/no_commit_maturity_2/poisson_9/poisson_d.pdf'), to = file.path(output_dir, 'appendix_plots/poisson_d_low.pdf'))
file.copy(from = file.path(output_dir, 'poisson_model/no_commit_maturity_2/poisson_5/poisson_g.pdf'), to = file.path(output_dir, 'appendix_plots/poisson_g_high.pdf'))
file.copy(from = file.path(output_dir, 'poisson_model/no_commit_maturity_2/poisson_9/poisson_g.pdf'), to = file.path(output_dir, 'appendix_plots/poisson_g_low.pdf'))

file.copy(from = file.path(output_dir, 'poisson_model/commit_maturity_2/poisson_13/debt_price_Poisson_model_commit.pdf'), to = file.path(output_dir, 'appendix_plots/debt_price_Poisson_model_commit.pdf'))
file.copy(from = file.path(output_dir, 'poisson_model/commit_maturity_2/poisson_13/value_function_Poisson_model_commit.pdf'), to = file.path(output_dir, 'appendix_plots/value_function_Poisson_model_commit.pdf'))

file.copy(from = file.path(output_dir, 'poisson_model/x_bar_vs_dt_transformed.pdf'), to = file.path(output_dir, 'appendix_plots/x_bar_vs_dt_transformed.pdf'))
file.copy(from = file.path(output_dir, 'poisson_model/target_x_vs_dt_transformed.pdf'), to = file.path(output_dir, 'appendix_plots/target_x_vs_dt_transformed.pdf'))

file.copy(from = file.path(output_dir, 'risk_aversion/2/cons_base_case.pdf'), to = file.path(output_dir, 'appendix_plots/risk_aversion_2_cons.pdf'))
file.copy(from = file.path(output_dir, 'risk_aversion/2/g_base_case.pdf'), to = file.path(output_dir, 'appendix_plots/risk_aversion_2_g.pdf'))
file.copy(from = file.path(output_dir, 'risk_aversion/2/d_base_case.pdf'), to = file.path(output_dir, 'appendix_plots/risk_aversion_2_d.pdf'))
file.copy(from = file.path(output_dir, 'risk_aversion/2/spread_base_case.pdf'), to = file.path(output_dir, 'appendix_plots/risk_aversion_2_spread.pdf'))
file.copy(from = file.path(output_dir, 'risk_aversion/risk_aversion_dist.pdf'), to = file.path(output_dir, 'appendix_plots/risk_aversion_dist.pdf'))
file.copy(from = file.path(output_dir, 'risk_aversion/risk_aversion_x_bar.pdf'), to = file.path(output_dir, 'appendix_plots/risk_aversion_x_bar.pdf'))
file.copy(from = file.path(output_dir, 'risk_aversion/risk_aversion_x_mean.pdf'), to = file.path(output_dir, 'appendix_plots/risk_aversion_x_mean.pdf'))
file.copy(from = file.path(output_dir, 'risk_aversion/risk_aversion_default_rate.pdf'), to = file.path(output_dir, 'appendix_plots/risk_aversion_default_rate.pdf'))

file.copy(from = file.path(output_dir, 'risk_aversion_vs_m/risk_aversion_sensi_m_cons_mult.pdf'), to = file.path(output_dir, 'appendix_plots/risk_aversion_sensi_m_cons_mult.pdf'))
file.copy(from = file.path(output_dir, 'risk_aversion_vs_m/risk_aversion_sensi_m_default_rate.pdf'), to = file.path(output_dir, 'appendix_plots/risk_aversion_sensi_m_default_rate.pdf'))
file.copy(from = file.path(output_dir, 'risk_aversion_vs_m/risk_aversion_sensi_m_x_bar.pdf'), to = file.path(output_dir, 'appendix_plots/risk_aversion_sensi_m_x_bar.pdf'))
file.copy(from = file.path(output_dir, 'risk_aversion_vs_m/risk_aversion_sensi_m_x_mean.pdf'), to = file.path(output_dir, 'appendix_plots/risk_aversion_sensi_m_x_mean.pdf'))

file.copy(from = file.path(output_dir, 'base_case_2/value.pdf'), to = file.path(output_dir, 'appendix_plots/value_2.pdf'))
file.copy(from = file.path(output_dir, 'base_case_2/g.pdf'), to = file.path(output_dir, 'appendix_plots/g_2.pdf'))
file.copy(from = file.path(output_dir, 'base_case_2/debt.pdf'), to = file.path(output_dir, 'appendix_plots/debt_2.pdf'))
file.copy(from = file.path(output_dir, 'base_case_2/spread.pdf'), to = file.path(output_dir, 'appendix_plots/spread_2.pdf'))